A Systematic Framework for the Development and Analysis of Signed

Aug 30, 2003 - The issue of systematic development of graph representations for chemical processes has not been addressed in the literature. This is a...
0 downloads 9 Views 510KB Size
Ind. Eng. Chem. Res. 2003, 42, 4789-4810

4789

A Systematic Framework for the Development and Analysis of Signed Digraphs for Chemical Processes. 1. Algorithms and Analysis Mano Ram Maurya,† Raghunathan Rengaswamy,‡ and Venkat Venkatasubramanian*,† Laboratory for Intelligent Process Systems, School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907, and Department of Chemical Engineering, Clarkson University, Potsdam, New York 13699-5705

In the recent past, graph-based approaches have been proposed by various researchers for safety analysis and fault diagnosis of chemical process systems. Though these approaches have shown promise, there are a number of important issues that have not been adequately addressed in the literature. The issue of systematic development of graph representations for chemical processes has not been addressed in the literature. This is an important issue because the development of digraphs is error-prone and time-consuming. Further, little attention has been paid toward understanding the conceptual relationship between the underlying mathematical description and the analysis procedures for the graph model. Also, the utility of these graphbased approaches at a flowsheet level has not been studied. With these issues in perspective, in this first part of the two-part paper, we focus on the systematic development of graph models and the conceptual relationship between the analysis of graph models and the underlying mathematical description of the process. The elimination of spurious solutions through the use of noncausal redundant equations and thorough analysis of inverse and compensatory response is also discussed. 1. Introduction Causal models are representations of cause and effect behavior in a system and hence are important in problems such as process fault diagnosis, hazard and operability analysis, operator training, and so on.1-5 Causal models are used in these applications to reason about the behavior of the process under various normal and abnormal conditions. Qualitative causal models in the form of graphs [directed graphs (DG) and signed DGs (SDG)] are often used. While the DG model captures the information flow, the SDG model captures both the information flow and the direction of effect (increase and decrease). Venkatasubramanian et al.6 have presented a comprehensive review of SDG-based methods and their applications in chemical engineering. To motivate the necessity of the research work presented in this paper, a succinct review is presented below. Iri et al.7 were the first to use SDG for fault diagnosis. In the following years, Umeda et al.8 showed how SDG can be obtained from differential algebraic equations (DAEs) for the process, though their methodology is limited to DAE systems having only one perfect matching so that the variable to algebraic equation (AE) matching is straightforward. An improved algorithm is used for applying the technique to a flowsheet level problem.9,10 Another group of researchers used partial system dynamics, statistical information about equipment failure, and digraphs to represent the failure propagation network for identifying fault location.11,12 * To whom correspondence should be addressed. Tel.: (765) 494-0734. Fax: (765) 494-0805. E-mail: [email protected]. † Purdue University. ‡ Clarkson University.

No sign is required in this analysis, but the method is limited to systems without any feedback. A rule-based method using SDG has been used for fault diagnosis by Kramer and Palowitch.13 In all of the above work, little attention has been given to the underlying relationships between the analysis of the mathematical model and the analysis of the digraph. Oyeleye and Kramer14 have presented important work on steadystate qualitative simulation. Algorithms for systematic use of steady state as well as dynamic model equations and process information have been given. The main focus of the work by Oyeleye and Kramer14 has been on eliminating spurious solutions, while ensuring completeness, through the use of noncausal and causal redundant equations. Though this work improved the understanding of systems exhibiting inverse response (IR) and compensatory response (CR), no formal proof has been presented as to why nodal balance on the extended SDG (ESDG) is valid. Chang and Yu15 have reported various techniques useful in simplifying SDGs for fault diagnosis. They show that ambiguities in control loops can be resolved by writing the controller equations in discrete and velocity form. Wilcox and Himmelblau3,4 have approached the problem of fault diagnosis using possible cause and effect graph models. The focus of their work is on efficient and complete fault diagnosis (by way of incorporating more information about the process, which could be quite specific sometimes) rather than promoting fundamental understanding of causal models. Vaidhyanathan and Venkatasubramanian16 have used digraph-based models for automated HAZOP analysis. HAZOP-digraph models are used to accomplish the task. SDG used in their work is generated from the steady-state description of the process as discussed by Mylaraswamy et al.17 (analysis of the validity of the algorithm has not been discussed).

10.1021/ie020644a CCC: $25.00 © 2003 American Chemical Society Published on Web 08/30/2003

4790 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

The recent work on SDG applications is by Rengaswamy and his colleagues.18-22 They use digraph models for optimal sensor location for efficient fault diagnosis. Palmer and Chung23 have used SDG-based models for HAZOP analysis. They have not made any attempt to check the validity of the SDG for individual units. Rather, they have tried to mitigate the problems arising as a result of using incorrect SDGs for individual units by proposing what is referred to as the modular approach in their equipment model builder. The most recent and promising application of SDG is the work presented by Chen and Howell24 on SDG-based analysis of control loops for fault diagnosis. On the basis of the above review, it is amply clear that, though the literature is replete with applications of graph models, the fundamental question of development and analysis of these models is not completely addressed. For example: (i) The process of developing these models is prone to human error and hence unreliable. There has been some effort toward developing these models systematically from the structure of the underlying mathematical description of the system.17 However, the correctness of the development strategy has not been addressed. (ii) The mathematical basis that underlies the various techniques that have been proposed in the literature for analyzing graph models has not been clearly articulated. (iii) Graph models would have a clear advantage over qualitative equations if one can use propagation through the graph to predict the behavior of the system. For example, fault diagnosis performed by propagating from the symptom nodes to the fault nodes using efficient graph theoretical search methods would be computationally much better than a generate-test method for solving qualitative equations. However, the validity of propagation to predict process behavior from a graph model has not been addressed adequately. It needs to be realized that, depending upon whether the SDG developed corresponds to a differential equation (DE), an algebraic, or a DAE system, the underlying properties could be completely different, and hence the same analysis methodology may not be valid for all SDGs. Thus, comprehensive, correct, and detailed analysis of various algorithms, methodologies, and issues in SDG analysis is the theme of this paper. It should be mentioned that lesser attention has been devoted to expedite the search process in fault diagnostic application, which actually has been the focus of many previous contributions in this area. The reason is that once a fundamentally consistent framework is developed, the existing search techniques could be applied, possibly after some modification. In this paper, development of graph models for various types of systems (DE and DAE systems and so on) will be discussed. The mathematical basis underlying the analysis of these models is then laid down. Clear emphasis is placed on distinguishing between transient and steady-state analysis. In this context, the idea of analyzing the graph models for the prediction of the initial response is introduced. Prediction of the initial response is valuable for incipient fault diagnosis and hence forms an important component of the analysis that is presented in this paper. The conditions under which propagation to predict process behavior would be valid are also discussed. Whenever appropriate, significant results will be supported through simple case studies from the literature. In the next section, the

Figure 1. Digraph for the DE system of example 1.

algorithms for digraph generation for different systems are presented. Section 3 deals with the understanding of the relationship between the algorithms and the corresponding analysis methodologies. Section 4 deals with the use of noncausal redundant equations for the elimination of spurious solutions. Section 5 deals with analysis of systems exhibiting IR and CR. 2. Algorithms for Digraph and SDG Generation The mathematical preliminaries for understanding the algorithms (dependency graph, bipartite graph, perfect matching, etc.) have been described in Appendix 1. In this section, algorithms for generating digraphs for systems described by DE, AE, and DAE are discussed. 2.1. Algorithm DE, for Systems Described by DEs. In systems described by DEs, explicit causality is from right to left.14 A one-step algorithm for the generation of DG is given below. For every DE, the variable on the left-hand side is matched with that equation and directed arcs are drawn from all of the variables on the right-hand side to the system variable on the left-hand side in that equation. Self-arcs also are drawn. If the expression for dxi/dt is not a function of xi, then a zero self-arc is drawn to xi. An illustrative example is presented below. Example 1. Let the system of AEs be

dx1/dt ) f1(e1,e2,x1,x3)

(1)

dx2/dt ) f2(e2,x1,x4)

(2)

dx3/dt ) f3(e1,e2,x1,x3,x4)

(3)

dx4/dt ) f4(e2,x2,x4)

(4)

The digraph for the above system is shown in Figure 1. 2.2. Algorithm AE, for Systems Described by AEs. Algebraic systems capture instantaneous behavior. Thus, there is no causality. Still, the set of system variables can be partitioned into several disjoint subsets [basically identifying strongly connected components (SCCs)] to make the information flow clear. Thus, causality is not explicit, and the digraph thus generated is noncausal in nature. The algorithm follows: [1] For every equation, directed arcs are drawn from exogenous variables to all system variables in that equation. Bidirected arcs are drawn among all system variables in that equation. The resulting digraph is called the dependency graph. [2] A bipartite graph between the equations and the system variables is drawn, and perfect matching is

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4791

Figure 2. Digraph for the AE system: various steps for digraph generation.

performed. Let the ith equation be matched with xj. Now all arcs from the exogenous variables in the ith equation to system variables other than xj are removed. Similarly, all arcs between xj and other system variables (xk*j’s) in this equation are made unidirectional toward xj, and all arcs among other system variables are removed. This is repeated for every equation. In step 2 of the above algorithm, any perfect matching can be selected because if more than one perfect matching exists, then there is at least one SCC and the SCCs consist of system variables that lead to multiple perfect matchings. Once the digraph is generated, each SCC can be replaced by a single node called a supernode to get a directed acyclic graph (DAG). In a DAG, there are no SCCs. The above procedure is illustrated for the steady-state equations derived from the DEs (by equating them to zero) in example 1. Figure 2a shows the dependency graph, Figure 2b shows the bipartite graph between the equations and the system variables, Figure 2c shows a particular perfect matching, Figure 2d shows the reduced digraph after deleting the arcs in step 2 of the algorithm, and Figure 2e shows the digraph after the SCC is replaced by a supernode. While using the above algorithm, if there is only one system variable in each equation, then information flow is from exogenous variables to the system variable present in the equation and there are no arcs among the system variables. Hence, there are no arcs that need to be deleted in step 2, and hence no solutions could be lost. When there is more than one system variable in some of the equations (a special case of this situation may be that more than one perfect matching and hence SCCs exist), step 2 of the above algorithm apparently means that only a particular system variable (which is matched with the equation) is affected by exogenous variables and other system variables in that equation. However, this is not true when SCCs exist. Hence, it has to be shown that none of the actual solutions are

lost while dropping some of the directed arcs in step 2 of the algorithm. Correctness of the algorithm has been shown in section 3. 2.3. Algorithm DAE, for Systems Described by DAE. The previous two algorithms (DE and AE) can be combined to develop the algorithm for DAE systems. Let the DAE system be given by

dX1/dt ) F1(xi∈X1,xj∈X2,el∈E) F2(xi∈X1,xj∈X2,el∈E) ) 0 where F1 and F2 are function vectors, X ) X1 ∪ X2, and X1 ∧ X2 ) φ. Thus, the set X is partitioned into two disjoint sets. The algorithm is as follows: [1] Apply algorithm DE on DEs and generate the first part of the digraph. [2] Apply algorithm AE on AEs, treating ∀xi ∈ X1 as pseudo exogenous variables, and thus generate the second part of the digraph. [3] Superpose the two parts to get the digraph for the entire system. The essential idea is to match the variables of the set X1 with DEs and those of the set X2 with AEs. As an illustrative example, consider the system consisting of eqs 1 and 2 and the AEs obtained from eqs 3 and 4 (by equating them to zero). It is easy to see that only one perfect matching exists between the set X2 ) {x3, x4} and the AEs (3) and (4). The resulting digraph is shown in Figure 3. The digraph thus generated corresponds to a dynamic model and hence can be used for predicting the initial response. To generate the digraph only for steady-state analysis, one can convert the DAE system into the AE system and algorithm AE can be applied. 2.4. Developing SDG from DG. In developing a DG, one need not consider the nonlinearity of the model because digraphs capture only information flow. On the other hand, SDG is used to predict qualitative departure from a given steady state. Thus, SDG also contains

4792 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Figure 3. Digraph for the DAE system.

information about the nature of the linearized system (around a given steady state), and the qualitative behavior is predicted correctly only in the qualitative linear region (i.e., the region in which the first nonzero partial derivatives do not change their sign) around that steady state. To achieve this, the given system is often linearized about a given steady state. Now all of the variables refer to differential values (i.e., X - X0, where X0 is a steady-state value) rather than absolute values, and they are called deviation variables. SDG can be obtained from DG by assigning appropriate signs to all of the directed arcs. Depending upon the nature of the model equations used to develop the DG (e.g., DE, AE, or DAE), arc signs are calculated as outlined below. 1. SDG for DE Systems. The sign of the arc from a variable on the right-hand side to the variable on the left-hand side is equal to the partial derivative of the expression on the right-hand side with respect to the chosen variable on the right-hand side. Mathematically, let the DE be

dxi/dt ) fi(xj∈X,el∈E) Then, the sign of the arc from xj to xi is

sign(xjfxi) ) [∂fi/∂xj] For the differential form of the equations to be valid for finite values of differential terms (e.g., dxj or del), the underlying assumption (to ensure conceptual consistency) is that ∂fi/∂xj is nonzero. If it is zero, then as mentioned in earlier literature,14 it should be replaced by ∂mfi/∂xjm if m is odd and by ∂mfi/∂xjm [dxj] if m is even, where m is the order of the first nonzero partial derivative. Similar arguments hold true for sign(elfxi). 2. SDG for AE Systems. If the ith equation is matched with xj and the other variables appearing in this equation are denoted by xk*j and el, then the sign of the arc from xk to xj is given by

sign(xkfxj) ) -([∂fi/∂xk]/[∂fi/∂xj]) ≡ [∂xj/∂xk] evaluated through the ith equation. If any of the first partial derivatives becomes zero, then the above expressions should be modified appropriately (similar to the SDG for DE systems). 3. SDG for DAE Systems. Recall that the digraph for DAE systems is generated by superposing the parts of the digraph generated using the algorithms DE and AE. Hence, the sign of the arcs in the resulting digraph

is also fixed using the rules outlined for DE and AE systems to generate the SDG. Once a SDG is generated, the next step is to design suitable method(s) to analyze the SDG to predict qualitative behavior of system variables for any deviation in exogenous variables (qualitative simulation) and vice versa (fault diagnosis). The next section deals with this issue. 3. Methods for SDG Analysis and Their Correctness When exogenous variables deviate from their nominal values, the system variables also deviate from their nominal values. The SDG developed in the previous section can be used to predict the sign of nodes representing system variables. As mentioned in the Introduction itself, correctness of the analysis methodology (in conjunction with the algorithm used for developing the SDG) is the most important issue. Given a SDG, regardless of the algorithm used for generating the SDG, intuitively one thinks that the sign of a node is governed by the total effect through arcs incident to that node (also known as nodal balance). Another way to predict the sign of a system variable could be to consider propagation through all of the directed paths from nonzero exogenous variables to the system variable. As shown later in this section, either of these might be wrong. Hence, the correctness of SDG solutions must be established. As we proceed, a number of definitions would be presented. First, we define nodal balance. Consider a SDG and a node in it representing a system variable xj. Let el and xk denote the exogenous variables and the system variables, respectively, from which there are directed arcs to xj. Definition 1. Nodal Balance. The qualitative state of a system variable is the net effect of all of the directed arcs incident on that node, i.e., the qualitative sum of the products of the signs of the incident arcs and the signs of the corresponding initial nodes. Mathematically,

[dxj] )

sign(xkfxj) [dxk] + ∑sign(elfxj) [del] ∑ k*j l

(5)

From the above equation, it is clear that one needs to know how to add and multiply two qualitative values. A very simple example is that the sign of the terminal node is the product of the sign of the initial node and the sign of the arc between them. As a consequence of this, if the sign of the arc from node xk to xj is + (-), then xj changes in the same (opposite) direction as xk. Rules for these qualitative operations (qualitative al-

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4793 Table 1. Laws of Qualitative Algebra series no.

[x]

[y]

[x] + [y]

[x][y]

δ([x])

1 2 3

0 ( (

[y] [x] -[x]

[y] [x] ?

0 + -

+ 0 0

gebra) have been discussed in the literature.7,14 Just for completeness, they are also listed in Table 1. In Table 1, “?” means undecided (i.e., the variables can have any sign) and δ([‚]) is the qualitative Dirac delta function. Next, the SDGs generated using the algorithms presented in section 2 are analyzed. Various analysis methodologies and the basis for their applicability have been discussed. 3.1. Analyzing SDG for DE Systems: Initial Response. SDG for the DE system is generated using the algorithm DE presented in the previous section. In this section, SDG is used to predict only the initial response, although useful causal information can be derived from this SDG to eliminate some of the spurious solutions in predicting the steady-state response of a dynamic system (section 5). First, the initial response is defined in a rigorous way. Then, we prove that propagation in SDG through the shortest paths predicts the correct initial response for DE systems. 3.1.1. Notion of the Initial Response. Given a dynamic system at a steady state, as one or more of the exogenous variables deviate from their normal operating values, the time derivatives of the system variables become nonzero and hence system variables also start deviating from their given steady-state values. Definition 2. Initial Response. The initial response of a system variable is its first nonzero response. Similarly, the initial response of the entire system is its response at the smallest time by which all of the system variables have shown their initial response. Definition 3. Shortest Path(s). Shortest path(s) from el to xj is (are) the directed path(s) with least length among all of the directed paths from el to xj. The above definition of the initial response would be compared with the traditional definition of the initial response presented in the literature.7,14 First we make the following observations. Observation 1. Relationship between the SDG and the Structure of the Matrices Describing a DE System. Let the DE system be described by

dX/dt ) F(xi∈X,el∈E)

(6)

Linearizing around the given steady state and using deviation variables, we get

dX/dt ) AX + BE

(7)

where A ) [aij] and B ) [bil] are the Jacobian matrices of the function F ) [fi] with respect to X ) [xj] and E ) [el], respectively, calculated as aij ) (∂fi/∂xj) and bil ) (∂fi/∂el). In the above expressions, X and E refer to the sets of deviation variables. Hereafter, the variables refer to deviation variables until otherwise mentioned. The context also should help in understanding the notation. It is worth mentioning that only step changes in E are considered but all of the analysis holds true for any kind of change (e.g., ramp, etc.) in E provided the sign does not change during the analysis. If a variable appears on the right-hand side in an equation and the corresponding partial derivative is

nonzero (if zero, the above expressions should be modified as suggested in section 2.4; for simplicity, we assume nonzero), then the appropriate entry in A, aij, is nonzero. A simultaneous look at the SDG shows that there is the directed arc (xj f xi) with the sign [aij]. We also know that the (i, j)th entry of an adjacency matrix for a given digraph is nonzero if there is a directed arc (directed path of length 1) from node i to node j in the digraph. Thus, the matrix A is essentially the transpose of the signed adjacency matrix of a part of the SDG containing all of the system variable nodes and the arcs among them. Similarly, B is the transpose of the signed adjacency matrix of a part of the SDG containing all of the arcs from exogenous variables to the system variables. Observation 2. Qualitative Equivalence of Propagation in a SDG and Matrix Multiplication of the Corresponding Matrices. Consider the matrices A and B with their elements in the symbolic form (zero entries are retained as 0 only) and the SDG for the DE system in eq 7. We analyze the expressions BE, ABE, A2BE, etc. The qualitative matrices [A] and [B] represent the directed arcs (the directed paths of length 1). Thus, the term [B][E], i.e., the qualitative equivalent of BE, represents propagation from the exogenous variables to the system variables through the directed path(s) of length 1. It is well-known that, given a digraph and its adjacency matrix A (a square matrix), if k is the smallest positive integer such that the (i, j)th element of Ak is nonzero (in the sense of symbolic value), then there is at least one directed path of path length k from the node xi to the node xj in the digraph. If self-cycles (if any) are not considered (i.e., ∀ i, aii ) 0), then the (i, j)th element of Ak is nonzero if there is at least one directed path of path length k from the node xi to the node xj in the digraph.25 Using the above results and recalling that the adjacency matrix for the SDG of the DE system in eq 7 is AT, one can verify for the DE system that if k is the smallest positive integer such that the (i, l)th element of Ak-1B is nonzero (in the sense of symbolic value), then there is at least one directed path of path length k from the node el to the node xi in the SDG. Thus, the qualitative equivalent of Ak-1BE, i.e., [A][A]...[A][B][E], represents propagation through the directed path(s) of length k from the exogenous variables to the system variables. Once matrices A and B are known, all of the directed paths from one node to another can be identified using any appropriate graph search method [such as depthfirst search (DFS)]. One can assign length of 1 unit to the directed arcs, and the shortest paths from exogenous variables to system variables can be identified to predict the initial response. Observation 3. Comparison with Traditional Definitions of the Initial Response. Traditionally, the initial response is defined to be the overall effect of propagation through all of the directed paths (also known as simple causal paths) from an exogenous variable to a system variable.14,26 Clearly, our definition of the initial response is more strict and would generate lesser spurious solutions for the initial response. We make the following claim. Claim 1. The initial response of a system variable xj due to changes in an exogenous variable el, predicted by propagation through all of the shortest path(s) from el to xj in the SDG for the DE system, is correct.

4794 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 Table 2. Propagation Table for the DEE (Example 2) exogenous variable system variable

S1

B1

B2

F

CF

hF

+ +

+ + +

W1 C1 h1 W2 C2

Results Obtained by Using Claim 1 0 + + 0 0 + 0 0 + + + 0 -

W1 C1 h1 W2 C2

Results Obtained by Using All of the Paths 0 ? ? + 0 0 + + 0 0 ? ? ? ? ? + 0 ? ?

+ + ? +

The proof of the above claim is presented next. 3.1.2. Correctness and Sufficiency of Propagation through Shortest Paths. The correctness of predicting the initial response can be established by analyzing the relationship between the qualitative equivalent of numerical dynamic simulation and propagation through the SDG. Consider qualitative dynamic simulation of the system described by eq 7. With X(0) ) 0, because all of the variables are assumed to be deviation variables, and step input in the exogenous variables, we have

X(t) ) eAtA-1(I - e-At)BE which simplifies to At2 A2t3 + + ... BE X(t) ) It + 2! 3!

(

)

}

(8)

Clearly, for a small enough time, the sign of X(t) will be the same as the sign of the first term in the series, and this will be the initial response of the system. On the basis of observation 2, one can see that a system variable would potentially show a nonzero sign if there are one or more directed arc(s) from one or more nonzero exogenous variables. Further, if there is no direct path (directed arc) from the exogenous variable to a particular system variable, then the initial response will be determined by the second term in the series, which is equivalent to directed paths of length 2 and so on. Hence, by following this argument, one can conclude that propagation through the shortest paths predicts the correct initial response. Thus, if the length of shortest path(s) from el to xi is k e n, then the initial response of xi is the first nonzero kth term in the series. For completeness of the proof, consider some special issues. Self-Cycles. As discussed above, in matrix multiplication, the effect can propagate through self-cycles only after reaching the node under consideration. Thus, selfcycles do not carry any information for predicting the initial response. For this reason (though not clearly mentioned), one can often see in the literature that selfarcs are not included in the SDG except for predicting IR and CR.14 IR and CR are considered in detail in section 5 in the context of reducing the number of spurious solutions in steady-state analysis of dynamic systems. Cycles. Propagation through cycles carries the feedback effect. However, this propagation starts only after the concerned system variable has shown a nonzero response. Thus, propagation through cycles can be used for predicting the transient behavior after the initial response. Example 2. Initial Response of a Double-Effect Evaporator (DEE). The mathematical model (a DE system) of

the DEE has been presented in the earlier literature.27 In the DE given below, the matrices A and B are shown as signed matrices only.

[ ] [ ][ ] [ ][ ] W ˙1 0 C˙ 1 0 h˙ 1 ) 0 W ˙2 0 C˙ 2 0 +

+ +

0 0 0 0 0

0 0 0 + -

W1 C1 h1 + W2 C2 0 0 + 0 0

0 0 + -

0 0 0 0

+ 0 0

0 + 0 0 0

0 0 + 0 0

S1 B1 B2 F CF hF

(9)

Though not shown here, the SDG for the above DE system can be easily drawn using the algorithm DE. Using claim 1, the initial response of all of the system variables for positive deviations in exogenous variables (single faults) is calculated, and the results are listed in Table 2. Such a table is called the propagation table. The results are interpreted as follows. When S1 becomes “high”, W1 shows a negative deviation during the initial evolution. One can see that the initial response for all of the system variables is unambiguous. As a comparison, the initial response is predicted by using propagation through all of the directed paths (traditional definition of the initial response), and the results thus obtained also are shown in Table 2. Table 2 shows that the signs of a number of system variables are ambiguous when all of the directed paths are considered. Thus, our definition of the initial response is more rigorous and better for SDG-based applications. 3.2. Analyzing SDG for AE Systems. Analyzing the evolution of a dynamic system involves three phases, namely, initial response, intermediate transients, and final response. The initial response has already been presented in reasonable detail. Analyzing intermediate transients involves propagation through cycles in the SDG, which does not provide any insight. The final response involves analyzing the corresponding AEs and is the subject matter of this section. The algorithm to generate SDG for systems described by AE (algorithm AE) has been presented in section 2.2. Now we prove the correctness of the algorithm AE and then discuss the equivalence of SDG and qualitative equations with respect to nodal balance followed by a discussion on the noncausal nature of the SDG for AE systems. Next, we present the results, which show that propagation through SDG is valid only when there is only one perfect matching and that in the case of multiple perfect matchings propagation through the SDG corresponding to only one of the perfect matchings is incomplete. We can make the following claim. Claim 2. The algorithm AE is correct, and SDG and qualitative equations are equivalent with respect to nodal balance. Proof of the above claim is discussed in sections 3.2.1 and 3.2.2. 3.2.1. Correctness of the Algorithm AE. The basic idea in establishing the correctness of the algorithm is to analyze various steps of the algorithm AE and then show that all of the steps are motivated by quantitative analysis. The validity of the algorithm is shown through the following steps:

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4795

1. The system of AEs is written in differential form. 2. Each of these equations is matched with exactly one system variable to get a perfect matching. (For any solvable AE system, at least one perfect matching exists.28) 3. Each equation can be rewritten in a slightly different form. In each equation, the variable that is matched with that equation is kept on the left-hand side and rest of the variables are kept on the right-hand side of the equation. 4. Now this representation can be converted into graphical form. Corresponding to a particular equation, directed arcs are drawn to the “node representing the variable matched with that equation” from the nodes representing other variables in that equation. Thus, solving the digraph is equivalent to solving the qualitative equations derived from a quantitative representation of the equations. The above steps have been explained below in the context of the algebraic system derived from the DE system of example 1. The equations in their differential forms are given below.

∂f1 ∂f1 ∂f1 ∂f1 de1 + de2 + dx1 + dx ) 0 ∂e1 ∂e2 ∂x1 ∂x3 3

(10)

∂f2 ∂f2 ∂f2 de2 + dx1 + dx ) 0 ∂e2 ∂x1 ∂x4 4

(11)

∂f3 ∂f3 ∂f3 ∂f3 ∂f3 de + de + dx + dx + dx ) 0 ∂e1 1 ∂e2 2 ∂x1 1 ∂x3 3 ∂x4 4 (12) ∂f4 ∂f4 ∂f4 de2 + dx2 + dx ) 0 ∂e2 ∂x2 ∂x4 4

(13)

On the basis of the perfect matching, eqs 10-13 are rewritten as follows.

∂f1 ∂f1 ∂f1 ∂f1 dx1 ) de1 de2 dx ∂x1 ∂e1 ∂e2 ∂x3 3

(14)

∂f2 ∂f2 ∂f2 dx4 ) de2 dx ∂x4 ∂e2 ∂x1 1

(15)

∂f3 ∂f3 ∂f3 ∂f3 ∂f3 dx ) de de dx dx (16) ∂x3 3 ∂e1 1 ∂e2 2 ∂x1 1 ∂x4 4 ∂f4 ∂f4 ∂f4 dx2 ) de2 dx ∂x2 ∂e2 ∂x4 4

(17)

Now a digraph can be generated, as discussed in step 4 of the procedure given above. Thus, we see that steps 2-4 together correspond to step 2 of the algorithm AE. In fact, there is no need to draw the dependency graph. Step 2 of the algorithm is sufficient for digraph generation. 3.2.2. Structural Equivalence of Qualitative Equations and the Digraph. The procedure for deriving the qualitative equations from steady-state quantitative equations is the same as that outlined by Oyeleye and Kramer.14 Qualitative equations can be obtained from the differential form of quantitative equations (eqs 10-13) by replacing the partial deriva-

tives and differential terms with their respective signs. The corresponding qualitative equations are given below.

[ ]

[ ]

[ ]

[ ]

∂f1 ∂f1 ∂f1 ∂f1 [de1] + [de2] + [dx1] + [dx3] ) 0 ∂e1 ∂e2 ∂x1 ∂x3 (18)

[ ]

[ ]

[ ]

∂f2 ∂f2 ∂f2 [de2] + [dx1] + [dx4] ) 0 ∂e2 ∂x1 ∂x4

[ ]

[ ]

[ ] [ ] [ ] [ ] [ ]

(19)

∂f3 ∂f3 ∂f3 [de1] + [de2] + [dx1] + ∂e1 ∂e2 ∂x1 ∂f3 ∂f3 [dx3] + [dx4] ) 0 (20) ∂x3 ∂x4

[ ]

∂f4 ∂f4 ∂f4 [de2] + [dx2] + [dx4] ) 0 ∂e2 ∂x2 ∂x4

(21)

By comparison of the above procedure with eqs 18-21, it is clear that there is a one-to-one correspondence between qualitative equations and the digraph (or rather the SDG). The procedure for digraph generation from quantitative equations (or, for that matter, SDG from qualitative equations) has already been discussed. Inversely, given the SDG, qualitative equations can be generated by applying nodal balance to all of the nodes representing system variables in the SDG.14 Thus, SDG and qualitative equations are equivalent with respect to nodal balance. The next section deals with a discussion on the noncausal nature of SDG for AE systems. 3.2.3. Noncausal Nature of the SDG for AE Systems. The systems described by AEs represent an instantaneous system. Variables affect each other (as depicted by the corresponding SDG) simultaneously. There is no temporal order in which the variables affect each other. In other words, time does not play any role. Similarly, the qualitative equations derived from the SDG for AE systems by applying nodal balance (if only the SDG is given) are noncausal in nature. This is somewhat contrary to what the SDG is meant for and that which we often see in the literature. Some examples from the literature are worth examining. We consider two representative examples in the next paragraphs. Iri et al. never discuss SDG for an instantaneous system.7 The algorithm presented by them can be used to develop digraphs for systems described by DEs and DAEs with only one perfect matching in the AEs portion (or for a given equation-variable matching). The existence of a single perfect matching makes the precedence ordering in the AEs of the DAE system trivial, and there are no feedback paths in the algebraic portion itself. Thus, strictly speaking, the resulting SDG (as a whole) is not causal in nature. Nevertheless, it is quite reasonable to say that the DE portion is causal in nature because causality is not disturbed by the presence of AEs. A complete discussion on AE systems with only one perfect matching is presented in the next section. Another example from the literature is investigated below. A methodology for including stability criteria (to eliminate spurious solutions) in qualitative simulation for predicting the steady-state response of a dynamic (causal) system described by DE has been presented by

4796 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Oyeleye and Kramer.14 Confluences (qualitative equations) are derived by applying nodal balance on the nodes of ESDG. They have referred to them as causal confluences because the ESDG (or the SDG from which ESDG is generated) represented a dynamic model. 3.2.4. AE Systems with Only One Perfect Matching. Precedence ordering for AE systems having only one perfect matching is transparent. Assuming unique quantitative solvability of the underlying system, we can say that there is at least one equation that contains only one system variable. This equation can be solved first, and the value of the corresponding system variable becomes known. Next, there is at least one equation that contains only one unknown system variable, which can be solved immediately. In this way the entire system can be solved sequentially. As a fact, perfect matching gives the complete sequence for solving the set of equations numerically. The same property is used in qualitative simulation as well. In terms of SDG, there are no SCCs in the SDG of such a system. It is intuitive to think that one should be able to predict the response of the system just by propagation through the SDG. It is indeed correct as stated in the claim given below. Claim 3. The response of AE systems with only one perfect matching can be predicted just by propagation through the SDG. This is because precedence ordering for AE systems having only one perfect matching is transparent. There are no SCCs in the SDG of such a system. Proof. Correctness of Propagation through SDG for an AE System with Only One Perfect Matching. Consider the steady-state system derived from eq 7. We get the following AE:

AX + BE ) 0

(22)

Equation 22 is uniquely solvable and has only one perfect matching. Hence, using suitable row permutations on A and B and suitable column permutations (possibly different from the row permutations) on A only, A can be converted into a lower triangular matrix with nonzero diagonal entries. The order of variables in X changes accordingly to ensure that the AEs themselves do not change as a result of column permutations. Further, when the resulting equations are divided by the respective diagonal entries, A becomes a lower triangular matrix, with all of the diagonal entries being 1. Also, the variables can be renamed so that the unique perfect matching is monotonic (i.e., the ith equation is matched with the new xi). Without loss of generality, one can assume that A already has the above structure, and eq 22 can be written as

X ) (I - A)X - BE

(23)

Though eq 23 is written in implicit form, it can be solved easily by forward substitution. Also, note that all of the diagonal entries of (I - A) are zero, and upon expansion, if a variable xj appears on the left-hand side of the equation, then it does not appear on the right-hand side. Some steps of the forward substitution procedure are

In eq 24, Bi is the ith row of B. The conditions over the indices ml and ml+1 of the term amlml+1 (the argument of the ∏ function) are interpreted as given below:

m1 ) i, ml > ml+1, and mj+1 ) k and the summation term on the last line and hence the compact underlined expression represent nothing but all of the directed paths of length j + 1 from exogenous variables to xi (observation 2) because all possible permutations (of indices) of length j satisfying the above conditions are chosen. The term (-1)j appears because the sign of the arc from node xj to xi is [-aij], etc. Similar expressions (though in different context) have been presented by Oyeleye and Kramer (see eq 8 in the respective paper).29 Thus, once one realizes that propagation is valid in numerical solution, it is correct for qualitative simulation as well (observation 2). Another way to understand the validity of propagation is to realize that (i) the inverse of a nonsingular lower triangular matrix is also a lower triangular matrix having diagonal entries as the inverse of the diagonal entries in the original matrix (retains the sign) and (ii) the solution of eq 22 is also given by (-A-1BE). This completes the proof. The importance of the above result is that efficient graph theoretical search algorithms (that rely on propagation) can be used for fault diagnosis. A simple example is presented below to explain how a perfect matching is derived. Consider the following AE system:

[ ][ ] [ ] [ ] 1 -2 3 x1 1 0 3 -7 0 x2 + 0 e1 ) 0 0 4 0 x3 0 0

(25)

By permuting rows 1 and 3, then permuting columns 1 and 2, and finally dividing the rows by 4, 3, and 3, respectively, the following system (note the lower triangular matrix A) is obtained.

[

][ ] [ ] [ ]

x2 4/4 0 0 0 0 x1 + 0 e1 ) 0 -7/3 3/3 0 -2/3 1/3 3/3 x3 1/3 0

(26)

The above equation (which is in the standard form, with all diagonal entries being 1) can also be written in its matched form as given below (note the strictly lower triangular matrix I - A). One could rename the

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4797

Figure 4. SDGs for example 3.

variables [x2, x1, x3] as [y1, y2, y3] for representational simplification.

[] [

][ ] [ ]

x2 0 0 0 x2 0 x1 ) 7/3 0 0 x1 + 0 e1 x3 2/3 -1/3 0 x3 -1/3

(27)

A SDG can be easily drawn (not shown here) using the above equation. From the above equation, clearly, [x2] ) 0, [x1] ) 0, and [x3] ) -[e1]. One can verify that the same results are obtained by propagation from e1 to the respective nodes xi in the SDG. 3.2.5. AE Systems with Multiple Perfect Matching. In most of the cases dealing with steady-state analysis, one ends up with an AE system having multiple perfect matchings. For numerical simulation, one resorts to inverting the matrix (for linear systems) or Newton’s method (or any other method of approximation for nonlinear systems). The SDGs (corresponding to any of the perfect matchings) of such a system contain SCC(s). For any given perfect matching, the AE system can be written in the matched form (as in eqs 14-17). As discussed in section 3.2.4 (see eq 23), suitable row permutations on A and B and suitable column permutations on A yield a block lower triangular matrix and blocks along the main diagonal represent SCCs (somewhat similar to the reduction of DE systems into block lower triangular form29). Because of the presence of the blocks along the diagonal, it is impossible to write an explicit expression (eq 24) for xi and hence propagation cannot be used for analyzing the SDG of such an AE system. In this section, we show through an example that propagation through the SDG corresponding to only one of the perfect matchings is incomplete for such systems. In certain situations, the missed solution might correspond to the stable steady-state response of the dynamic system (described by DE) from which the underlying AEs were derived. We present such an example. Example 3. Consider the system given below.

}[

dx1 ) -7x1 - 7x2 + e1 dt dx2 ) x2 + 9x3 w dt dx3 ) 4x1 - 13x3 dt -7 - 7 0 1 9 and B ) 0 A) 0 1 4 0 -13 0

]

[]

(28)

The SDGs for the DE in eq 28 and the corresponding

AE are shown in Figure 4. The system described by eq 28 is a stable system. Further, as long as the qualitative nature of A is retained, for stable A there is a unique solution X ) [-, +, -]T for e1 ) +. Simply solving the qualitative equations (derived from AEs or derived by nodal balance on any of the SDGs for the AE) yields two solutions: the unique stable solution and an unstable solution. The unstable solution indeed corresponds to an unstable system that has its system matrix, say A′, from the same qualitative class as A (i.e., [a′ij] ) [aij], ∀ i, j). As an example, let a′ij ) aij, ∀ i, j, except that a′12 ) -2. Then A′ is unstable, and X ) -A′-1B ) [+, -, +]T is the unstable solution. Interestingly, the AEs corresponding to both the stable system and the unstable system result in the same qualitative equation. Thus, given the level of information, truly speaking, none of the solutions generated by the AEs are spurious. If more numerical information is available, one can derive some redundant equations (this could be based upon stability criteria as well) and use them to eliminate the unstable solution. Redundant equations are discussed in detail in sections 4 and 5. Now we analyze the SDG for the AE corresponding to the first perfect matching (see Figure 4b). Propagation through the SDG yields X ) [+, -, +]T for e1 ) + (the unstable solution). One can easily verify that propagation through this SDG never predicts the stable (correct) solution. However, propagation through the SDG for the AE corresponding to the second perfect matching does yield the stable solution (see Figure 4c). Still, considering only the AE, none of the SDGs (Figure 4b,c) yield a complete set of solutions individually by propagation. The reason is that, for AE systems with multiple perfect matchings, nodal balance and propagation through the SDG corresponding to only one of the perfect matchings are not equivalent. If one still wants to use propagation for SDGs for an AE system with multiple perfect matchings, then completeness can be ensured by propagating through all of the SDGs for the AE system corresponding to different perfect matchings. For example, all of the solutions predicted by propagation through SDGs of Figure 4b,c generate the complete set of solutions for the AE system. At the same time, as the AE system becomes more complex, the number of perfect matchings can grow combinatorially, and storing all such SDGs and analyzing them by propagation might not be better than the generate and test method for solving steadystate qualitative equations. The generate and test procedure works by generating all possible signs for variables and then testing to see which sign patterns satisfy the qualitative equations.29 Though the SDG cannot be used directly, its topological structure (e.g., the nodes that form a SCC) can be used to identify a partial order in which the signs of the variables should

4798 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Figure 5. SDGs for example 4.

be chosen and the order in which the equations should be tested to expedite the search process.29 Some more details on the generate and test procedure are provided in part 2 of this series of papers. Before we present another example, let us look at the SDG shown in Figure 4a. Propagation generates the correct initial response (as expected) but fails to predict the correct final response (SDGs for the AE contain a positive cycle). The dynamic system representing this SDG exhibits IR. Now we present an example of a reaction system. Example 4. An Autocatalytic Reaction System.14 The governing equations for the isothermal, irreversible, autocatalytic reaction A f B are given below:

dcA F(cA,i - cA) ) - krcAcB2 dt V

(29)

FcB dcB ) krcAcB2 dt V

(30)

In the above system, cA,i is the exogenous variable and cA and cB are the system variables. The AE system corresponding to the above DE system has two perfect matchings, and the two SDGs are shown in parts b and c of Figure 5, respectively. As discussed by Oyeleye and Kramer,14 cA exhibits IR with respect to cA,i (Figure 5a). One can easily verify that propagation through the SDG of Figure 5c yields the correct results (e.g., [cA, cB] ) [-, +] for cA,i ) +). One cannot get the actual result by propagating through the SDG of Figure 5b. This example clearly shows that one must be cautious in analyzing a SDG by using propagation if the SDG corresponds to an AE system with multiple perfect matchings. Essentially, for AE systems with multiple perfect matchings, the qualitative equations can be solved by using a generate and test procedure (see the exception below). There is one exception to the above result; viz., AE systems with multiple perfect matchings having only negative cycles in (any of) their SDGs. We make the following claim. Claim 4. If an AE system has multiple perfect matchings but the SDG (corresponding to a perfect matching) contains only negative cycles, then any perfect matching (for which the SDG contains only negative cycles) can be chosen and the response of the system can be predicted by propagation through all of the directed paths in the corresponding SDG. As an example, consider an AE system with two perfect matchings (and hence a negative feedback loop in the SDGs). One can see that propagation through the directed path(s) (not through the cycle; i.e., a negative feedback cycle can be teared) in any of the two SDGs

generates the same result. Thus, any of the SDGs can be chosen, and propagation yields the correct results. A proportional control loop in a static process is an excellent example of such a system. The proof of claim 4 is presented below. Proof. We revisit eqs 22 and 23. For an AE system with multiple perfect matchings, by suitable row permutations on A and B and suitable column permutations on A only and then by division of the resulting equations by the respective diagonal entries of A, A can be converted into a block lower triangular matrix, with all diagonal entries being 1 (B changes accordingly). One can realize that for AE systems with multiple perfect matchings the above transformation is not unique and eq 23 represents just one of the perfect matchings. In fact, once one such perfect matching is generated, other perfect matchings can be generated through suitable column permutations (which are related to the cycles and/or SCCs in the SDG) on A (standard form with ∀ i, aii ) 1) and by division of the rows by the corresponding diagonal entries (to get the standard form). In eq 23, (I - A) is also block lower triangular, its all diagonal entries are zero, and the negative sign before A and B accounts for the sign of the arcs in the corresponding SDG for the AE system (see section 2.4). The blocks along the main diagonal correspond to the SCCs in the SDG. Now we consider the imaginary dynamic system given below:

dX/dt ) -AX - BE

(31)

We make the following observations. Observation 4. Comparison of the SDG for the AE System and the Corresponding Imaginary DE System. The SDG for the DE system can be derived from the SDG for the AE system by adding a negative self-cycle at each node (because all of the diagonal elements of A are 1; such a DE system cannot exhibit CR14) representing a system variable. Further, because the SDG for the AE system contains only negative cycles, the SDG for the DE system also contains only negative cycles (and self-cycles). Observation 5. Ultimate Response of a Dynamic System. The ultimate response of the DE system represented by eq 31 is governed by eq 22, which represents the AE system of our interest. Further, it is a wellknown result that if the SDG of a DE system contains only negative cycles (including self-cycles), then its steady-state response can be predicted by propagation through all of the directed paths.29 It can be noted that Oyeleye and Kramer have used only a necessary condition for stability of a dynamic system (viz., (-1)n|A| > 0 for the DE system of eq 7),

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4799

which is satisfied by the DE system represented by eq 31.29 Also, if the SDG of a DE system contains only negative self-cycles (i.e., all of the diagonal elements of the system matrix are negative), then there exists a stable DE system with the same SDG (derived from theorem 1 in Basset et al).30 Since we are interested only in qualitative analysis, the imaginary DE system need not be stable for our analysis (observation 5) to be correct. By using observations 4 and 5, it is trivial to conclude that claim 4 is correct. This completes the proof. In the above claim and proof, one can note the stress on choosing a SDG with only negative cycles. Apparently, if the SDG for a particular perfect matching contains only negative cycles, then so do the SDGs for any other perfect matching. Another subtle observation is that different imaginary DE systems corresponding to eq 22 (or eq 23) in standard form for different perfect matchings need not have the same qualitative structure [i.e., not reducible to one another through row permutations and corresponding column permutations (i.e., relabeling the variables)]. Incompleteness of propagation through only one of the SDGs for an AE system with multiple perfect matchings (excluding the exception) prompts one to question the validity of propagation through only one of the SDGs for a DAE system with multiple perfect matchings in the AE part (part 2 of the SDG). Next DAE systems are analyzed. 3.3. Analyzing SDG for DAE Systems. Having analyzed DE and AE systems, DAE systems can be analyzed using the methods presented in sections 3.1 and 3.2 with minor modifications. First, we discuss the initial response for DAE systems with only one perfect matching. Then, we discuss issues involved in predicting the initial response of DAE systems with multiple perfect matchings. Finally, a brief discussion on the steady-state analysis of DAE systems is presented. The model of the DAE systems used for analysis is derived from the model used in section 2.3 by linearization.

dX1 ) A11X1 + A12X2 + B1E dt

(32)

A21X1 + A22X2 + B2E ) 0

(33)

3.3.1. DAE Systems with One Perfect Matching: Initial Response. The initial response of DAE systems with only one perfect matching can be predicted correctly by propagation through shortest paths by introducing some modifications in the definition of the path length in the SDG. To guide the proof and to identify the necessary modifications in the definition of the path length, eqs 32 and 33 are analyzed in the numerical domain. The analysis presented in this section is similar to the one presented in section 3.2.4. Assuming A22 to be a lower triangular matrix with diagonal entries 1 and hence nonsingular, eq 33 can be written in two different forms as shown below.

X2 ) -A22-1A21X1 - A22-1B2E

(34)

X2 ) (I - A22)X2 - A21X1 - B2E

(35)

and

For numerical treatment, eq 34 can be used for direct substitution for X2 in terms of X1 and E into eq 32 and

all of the methods discussed in section 3.1 can be used for the resulting DE system. However, in qualitative analysis, matrix inversion cannot be expressed in a straightforward way. Hence, eq 35 is used instead. This equation is very similar to eq 23. Thus, using the forward substitution, X2 can be written in terms of X1, E, A21, A22, and B2. As mentioned in section 2.3, xj ∈ X1 can be treated as a pseudo exogenous variable. Mathematically, one can write

(X2)i ) -(B2)iE i-j

j|mj+1)k

k)1

l)1|m1)i,ml>ml+1

[∑ (∑{[ ∏ [∑ (∑{[ i-1

(-1)

j)1

j

i-1

(A21)iX1 -

(-1)

j)1

] })]

i-j

j|mj+1)k

k)1

l)1|m1)i,ml>ml+1

j

E-

(A22)mlml+1 (B2)k



]

(A22)mlml+1

(A21)k

})]

X1 (36)

where (X2)i is the ith variable of the subset X2, (B2)k is the kth row of B2, etc. Various terms in eq 36 show the effect of directed paths from E and X1 to X2, and thus it can be used to obtain the solution of X2 qualitatively. Hence, propagation is correct. Next we show how propagation with a modified definition of shortest path(s) can be used for predicting the initial response for such a system. One can substitute for X2 (using eq 34) in eq 32, and the following expression is obtained.

dX1/dt ) [A11 + A12(-A22-1A21)]X1 + [B1 + A12(-A22-1B2)]E (37) The above system is similar to the system described by eq 7 with a change of notation

X ≡ X1, A ≡ A11 + A12(-A22-1A21), B ≡ B1 + A12(-A22-1B2) The matrices A and B can be seen as the sum of two matrices. Further, the matrices A12(-A22-1A21) and A12(-A22-1B2) are the expressions giving information about additional directed paths from the nodes xk ∈ X1 and el ∈ E to the nodes xi ∈ X1 through the nodes xj ∈ X2 (observation 2). Because these terms are added with A11 and B1, respectively, which give information about directed arcs (path length 1, observation 1), each term of these matrices too expresses a directed path with effective path length 1. Now different parts of the expression A12(-A22-1A21) are analyzed in the reverse order. A21 contains information about the arcs from xi ∈ X1 to xj ∈ X2 while constructing part 2 of the digraph in the algorithm DAE. The information related to the solution of X2 in terms of X1 and E (or propagation in the AE part of the SDG) is contained in the term -A22-1, and A12 contains information about the arcs from X2 to X1 while constructing part 1 of the digraph. Because of the causal nature of DEs, we say, A12 retains information about directed arcs (paths of length 1). Thus, if we assume that the arcs (or the directed paths) corresponding to A21 and -A22-1 have path length zero, then our statements about predicting the initial response (section 3.1) are still correct. Hence, we propose the following definition and claim.

4800 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Figure 6. Propagation for the DAE systems with one perfect matching.

Definition 4. Arc Length in SDG of a DAE System. The arcs in part 1 of the DG (algorithm DAE) have length 1. The arcs in part 2 of the DG have arc length 0. Claim 5. With the above definition, the initial response of a DAE system with only one perfect matching can be obtained using claim 1. The above claim has been already proved, and it is also explained in Figure 6. Now we present an example. Example 5. Exothermic Irreversible Continuous Stirred Tank Reactor (CSTR). Two reactants A and B are fed to a CSTR. The total volumetric flow rate is q. Inlet concentrations are cA,i and cB,i, respectively, and the inlet temperature is Ti. The second-order exothermic reaction A + B f P takes place inside the CSTR. The reactor jacket temperature is Tj. The mathematical model is given below.

dcA q(cA,i - cA) ) - k2cAcB dt V

Figure 7. SDG for the DAE system in example 5.

(38)

exogenous variable

dcB q(cB,i - cB) ) - k2cAcB dt V

(39)

dcP -qcP ) + k2cAcB dt V

(40)

dT q(Ti - T) (-∆H)k2cAcB hs(T - Tj) ) + dt V FcP VFcP k2 ) a2 exp(-E2/RT)

Table 3. Propagation Table for the CSTR (Example 5) system variable

(41) (42)

The SDG for the above DAE system (with one perfect matching) is shown in Figure 7. The exogenous variables are q, cA,i, cB,i, Ti, and Tj, respectively. The system variables are cA, cB, cP, T (T > Ti), and k2, respectively. The initial response of the system variables for changes in exogenous variables can be predicted by using claim 5. The results of propagation through the shortest paths are given in Table 3. One can see that the signs of all of the system variables (initial response) can be determined unambiguously. Next, the initial response is predicted by using all of the directed paths from the exogenous variables to the system variables. These

cA,i

cB,i

q

Ti

Tj

+ + +

+ + +

cA cB T cP k2

Results Obtained by Using Claim 5 + + + + + + + + + + -

cA cB T cP k2

Results Obtained by Using All of the Paths + ? ? + ? ? ? ? ? + ? ? ? ? ? ? ? +

? ? + ? +

results also are listed in Table 3. One can easily see that the results obtained by using all of the directed paths are very poor. 3.3.2. DAE Systems with Multiple Perfect Matchings. Incompleteness of propagation through the AE part of the SDG for a DAE system with multiple perfect matchings makes it prohibitive to use propagation through the entire SDG to predict the initial response of such a system. Two possible ways to overcome this problem are as follows: 1. Conversion of DAE into DE by Differentiation of the AEs. The AEs can be differentiated with respect

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4801

to time, and linear expressions for dX1/dt can be used to derive an expression for dX2/dt (treat X1 as a pseudo exogenous variable). This would require inverting the matrix A22, so this has to be carried out numerically. There is no nice way to carry out this calculation in the qualitative domain. Another problem might arise in differentiating the exogenous variables. Often one assumes a step change in exogenous variables. The time derivative of a step change requires modeling of the Dirac delta function in the qualitative domain. In the context of numerical simulation, this does not impose any problem because the variables receiving an arc from such exogenous variables show a finite change at time t ) 0+, which is not a serious problem. This cannot be handled in the qualitative domain. Still, one can predict the initial response qualitatively by assuming just a ramp (or first order) change in exogenous variables in the beginning. Once the system is in the form of a DE, all of the methods discussed previously for DE systems can be applied. One can note that the system matrix of the resulting DE system is singular because the later rows are linear combinations of the rows corresponding to dX1/dt. One should also note that the above procedure requires numerical treatment, and hence it renders SDG-based analysis useless except that the topological structure of the SDG can be used wherever beneficial. 2. Conversion of DAE into DE by Substitution for Algebraic Variables in the DEs. By substitution for X2 in terms of X1 and E (eq 34) in the DEs, the resulting system is again a DE system for X1. All of the methods discussed for DE systems are now applicable to the resulting system. The price is that substitution has to be carried out numerically. As suggested in section 3.2.5, an exception to the above result is the DAE systems with multiple perfect matchings containing only negative cycles in part 2 (corresponding to the AE portion) of their SDGs. Propagation can be used to predict the initial response of such a system. A proportional-integral control loop in a static process is a good example of such a DAE system. 3.3.3. Steady-State Analysis of DAE Systems. The steady-state response of a DAE system is governed by the corresponding AEs, and all of the discussion presented previously about AE systems can be used. In general, the number of perfect matching(s) does not decrease when a DAE system is converted into the corresponding AE system. For example, if the original DAE system has more than one perfect matching, then in most of the cases, the steady-state system also would have more than one perfect matching. To understand certain exceptions, the structural transformations that occur when a DAE system is converted into the corresponding AE system are investigated in the next paragraph. The perfect matching in the DE portion of the DAE system is trivial. However, when DEs are converted into AEs, the underlying perfect matching can change considerably. The most obvious and frequently encountered situation is the one in which one of the variables has a zero self-cycle (in the DE portion of SDG for the DAE system). The corresponding AE can no longer be matched with this variable. This, in turn, changes the equation-variable matching for other equations and variables as well. As studied by Oyeleye and Kramer,14,29 in such systems some other variable (other than the one having the self-cycle) can potentially show CR. In such cases, the number of perfect matching(s) in the

resulting AE system can decrease. In another situation, when there are SCCs in the DE portion as well as the AE portion of the SDG, the number of perfect matchings can increase to a great extent and the system can become even more complex for qualitative analysis as well as for numerical analysis (recall that SDGs capture the structural complexity of the system). In any case, the generate and test method can be used to predict the final response. 4. Reducing the Number of Spurious Solutions The prediction of spurious results is inherent in all qualitative reasoning based methods because some of the quantitative information is lost in the qualitative domain. Though the digraph solutions always include the actual solution, many spurious solutions are generated that reduce the utility of the digraph analysis. For example, in fault diagnosis, the set of root-cause nodes should be as small as possible. Thus, just the digraph is not very useful for qualitative analysis. Hence, there is a need for developing techniques that help in identifying and removing spurious solutions to make graphbased (such as SDG) analysis meaningful. In the context of qualitative dynamic simulation, Kuipers has shown that qualitative simulation can generate all possible solutions that can be exhibited by a real system and that the spurious solutions (generated as a result of the incompleteness of knowledge captured through the qualitative model) cannot be eliminated completely.31 Several other researchers have contributed to the work on elimination of spurious solutions in the context of qualitative dynamic simulation and semiquantitative or fuzzy-qualitative dynamic simulation.32-43 Because the focus of this paper is on the elimination of spurious solutions in SDG-based steady-state analysis, a formal and detailed discussion on these contributions is not presented. In section 3, it has been stressed that SDGbased analysis is mainly useful for the prediction of the initial response and steady-state response. Most of the contributions for the elimination of spurious solutions in SDG-based steady-state qualitative analysis can be attributed to Oyeleye and Kramer14,29 in which they use steady-state equations as well as the topological structure of the SDG (of the dynamic system) to generate redundant equations. By incorporation of more information, the number of spurious solutions in the prediction of the final response can be reduced. Additional information, broadly categorized as noncausal and causal information, can be used to generate noncausal and causal redundant equations, respectively. In this section, various methods for the generation of noncausal redundant equations have been presented. The generation and use of causal redundant equations is the focus of section 5. 4.1. Noncausal Redundant Equations. Noncausal redundant equations are generated from noncausal information, more specifically AEs describing the original steady-state system.14 Notice that this cannot be used for eliminating spurious solutions for predicting the initial response except in the AE portion of a DAE system. This section deals with generating such redundant equations at the unit model level. The generation of redundant equations at the flowsheet level has been discussed in Part 2 of this series of papers, where an overall framework for SDG-based flowsheet-level analysis has been presented. The usefulness of redundant equations for eliminating spurious solutions can be understood through the following example.

4802 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Figure 8. SDG analysis in example 6. Table 4. Results of Quantitative Analysis in Example 6 value of the variable description

e1

e2

x1

x2

x3

normal quantitative value new quantitative value new qualitative value

3.0 3.0 0

2.0 3.0 +

2.5 3.0 +

-0.7 -0.6 +

2.1 1.8 -

Table 5. Results of Qualitative Analysis for e2 ) + (Example 6) qualitative value description of information used

x1

x2

x3

without redundant equation with redundant equation

+ +

? +

-

Example 6. Consider the system

2x1 - e1 - e2 ) 0

(43)

x1 - 2x2 + x3 - 2e1 ) 0

(44)

2x2 + 4x3 - 3e1 + e2 ) 0

(45)

For the above example, SDG is shown in Figure 8. Table 4 shows the results of quantitative analysis for a change in the exogenous variable e2 from 0 to +. From the SDG of Figure 8, x1 ) +. It is easy to show that x3 ) because x3 ) + or 0 means x2 ) +, which in turn means x3 ) -, contradicting the assumption. To be precise, whenever two or more arcs with conflicting effects are incident on a node, all possible signs (+, 0, and -) are chosen for that node and tested against the corresponding qualitative equations derived by applying nodal balance. The signs of all of the nodes in a SCC are decided simultaneously so that all possible combinations of their signs are chosen and tested. In this way, one finds that x2 can have any qualitative state. If we use the redundant equation (46) (obtained by the manipulation [eq 43] - 2[eq 44] + [eq 45]), it gives one more qualitative equation, which sets x2 ) +, thus eliminating all spurious solutions. These results are summarized in Table 5.

6x2 + 2x3 ) 0

(46)

4.2. Methods for the Generation of Noncausal Redundant Equations. As discussed in the previous section, the basic idea behind the generation of noncausal redundant equations is the use of additional information and the original AEs. Still, these methods can be categorized based upon the types of equations and the type of additional information used. Various methods of generating noncausal redundant equations are discussed below.

4.2.1. Algebraic Manipulation. Redundant equations can be derived through algebraic manipulation of two or more equations of the original system.14 A redundant equation helps in reducing the number of spurious solutions if it contains a subset of variables different than the subset of variables present in the equations from which the redundant equation is derived. Consider a set of two equations having at least one system variable in common. Let the number of system variables in two equations be n1 and n2, respectively, and k (min(n1, n2) g k g 1) system variables be common in both. Then the total number of distinct system variables in this system of two equations is n1 + n2 - k. Using these two equations, one variable can be eliminated at a time. Hence, k redundant equations, each having a different subset of system variables, can be obtained, and the number of system variables in each redundant equation is n1 + n2 - k - 1 g max(n1, n2) 1. The new redundant equation contains a different subset of system variables, and hence it can be helpful in reducing the number of spurious solutions.14 4.2.2. Use of Numerical Information. Numerical information related to the physical and chemical properties of the species involved in the system can be used to reduce the number of spurious solutions. Such information can be included in the physical and chemical descriptions of the system. Consider the expression for the specific enthalpy of a liquid mixture of two components. The specific enthalpy increases with increases in the temperature and mole fractions of the more volatile component and less volatile component, respectively. Using the condition that mole fractions of two species add up to 1 and using the numerical values of the parameters involved in the corresponding quantitative equation, one can get a qualitative equation that has only two variables, viz., temperature and mole fraction of the more volatile component. Such an equation can be used in reducing the number of spurious solutions. 4.2.3. Expressions in Terms of Intensive Variables. Some of the system variables can be expressed in terms of exogenous variables and intensive system variables (e.g., mole fraction). This is another form of the redundant equation and can be included explicitly in the description of the unit model. As a simple example, consider the two equations given below describing molar balance in a flash vaporizer.

f)w+v fxf,1 ) wxw,1 + vxd,1 From these two equations, simple algebraic manipulation yields the expressions given below.

xd,1 - xf,1 w)f xd,1 - xw,1 xf,1 - xw,1 v)f xd,1 - xw,1 4.2.4. Use of the Relative Order of Magnitude. The relative magnitude of deviations in different system variables can be used in reducing spurious solutions. However, such equations cannot be stored in the model library because these are quite specific and hence the generality of the approach might be lost sometimes.

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4803

Another way to use the order of magnitude reasoning is by attaching a weight besides the sign to the directed arcs. Bhushan and Rengaswamy18 have used such an analysis for designing sensor networks to improve fault diagnostic resolution. 5. Generation of Causal Redundant Equations The generation of causal redundant equations depends on stability conditions that are imposed on the system matrix corresponding to the DEs. DEs model the temporal behavior (causal in nature) of a system, and hence the redundant equations thus obtained are also causal in nature. Oyeleye and Kramer14,29 introduced the idea of causal redundant equations and have discussed the generation of these equations in considerable detail. Causal redundant equations provide a powerful way in which spurious solutions could be removed, and hence we focus on this technique in this paper. In this section, we discuss some modifications/ corrections to the work presented by Oyeleye and Kramer.14 First, we explain the concept of causal redundant equations through a simple example and summarize the work done by Oyeleye and Kramer.14 For a detailed treatment and analysis of causal redundant equations, an interested reader is referred to work by Oyeleye and Kramer.14,29 Example 7. Consider the dynamic system dX/dt ) AX + BE and the corresponding steady-state system AX + BE ) 0 with the matrices A and B as given below.

[

]

[]

-7 -7 0 1 A) 0 1 9 and B ) 0 4 0 -13 0

(47)

We solve the steady-state system by algebraic manipulation. Consider the symbolic form of the matrices A and B in eq 47. Perform the following manipulation on the augmented matrix [A|B].

row(2) r row(2) + row(1) (-a22/a12) row(2) r row(2) + row(3) (-[row(2)]1/a31) where [row(2)]1 is the first element of the second row. The resulting equation is given below.

(

a23 +

) (

)

a33a11a22 a22 x3 + b )0 a31a12 a12 1

The sign of the coefficient of x3 can be calculated by using stability conditions. The stability of the dynamic system described by eq 47 requires that all of the eigenvalues of A have a negative real part. A necessary condition for the same is that all of the coefficients of the characteristic equation of A, viz., |λI - A| ) 0, must be positive and hence

(-1)n|A| > 0 By noting the signs of elements of A and B and that 0 > |A| ) a11a22a33 + a23a31a12 for stability, one can immediately conclude (using other qualitative equations as well and then solving) that [X] ) [-, +, -]T. This is the stable solution. In the above procedure, one can get some insight about what manipulations to perform. For example, in the SDG for the above system (not shown here), there

is a positive self-cycle on x2 and it is the second row on which the above elementary matrix operations are performed. As Oyeleye and Kramer14,29 have shown, a well-defined structure for identifying such manipulations is possible through the identification of IR and CR and hence the inverse variables (IVs) and compensatory variables (CVs) in the SDG. Concise definitions and descriptions of IR, IV, CR, and CV (as described by Oyeleye and Kramer14) are presented below. A process variable is said to exhibit IR with respect to an exogenous variable (and hence it is called an IV) if its final response is opposite to its initial response for nonzero deviations in the exogenous variable. Similarly, a process variable is said to exhibit CR with respect to an exogenous variable (and hence it is called a CV) if it shows zero steady-state deviation for nonzero deviations in the exogenous variable. IR and CR are the net result of opposite causal effects that arise from (i) multiple acyclic (feed-forward) directed paths to a system variable and (or) (ii) negative feedback loops. Thus, propagation through all of the (feed-forward) directed paths accounts for IR and CR arising from step i but fails to account for IR and CR arising from step ii. 5.1. Previous Work on Causal Redundant Equations. The results of the work done by Oyeleye and Kramer14 are briefly summarized below. (i) IR and CR in Dynamic Systems. The steady-state response of a stable dynamic system represented by DE can be analyzed by decomposing the system into SCCs and then solving the variables of each SCC with respect to their local disturbances. A necessary condition for IV and CV (with respect to a local exogenous variable) is that (1) IV/CV should be located inside a SCC on a negative feedback cycle. (ii) Necessary Conditions for IR. Other necessary conditions for IR are as follows: (2) The complementary subsystem to at least one of the acyclic paths from the local exogenous variable to the IV should contain a positive cycle (or self-cycle). (3) The complementary subsystem to the positive cycle in one of the complimentary subsystems in condition 2 should not have zero signed determinant (i.e., violating conditions 4a, 4b, or 5 discussed below). (iii) Necessary and Sufficient Conditions for CR. The other necessary and sufficient conditions for CR are as follows: (4) The complementary subsystems to all acyclic paths from the local exogenous variable to the CV should each (a) have at least one zero self-cycle and (b) not have a cycle containing all of the variables in the subsystem. (5) The complementary subsystems to all nonzero cycles (excluding self-cycles) in each complimentary subsystem in condition 4 should each satisfy conditions 4a, 4b, and 5. (iv) Generating Causal Redundant Equations from ESDG. After identification of IV/CV, ESDG is drawn by adding additional arcs in the SDG for the DE system appropriately. IVs are located as strings on the negative feedback cycles. For each string of IVs, a directed arc is drawn from the variable located just upstream of the first IV in the string to the noninverse variable located just downstream of the last IV in the string. The sign of this arc is the same as the sign of the directed acyclic path from the upstream variable to the downstream variable. Similarly, for each string of CVs, directed arcs are drawn. The only difference is that the sign of the arc is (δ([dxj]), where [dxj] is the sign of last CV on the CV string. Once ESDG is drawn, nodal balance is

4804 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

applied on the nodes of the ESDG to get the desired causal redundant equations. Having discussed the existing work on IV/CV by Oyeleye and Kramer,14,29 a number of modifications to their work are presented. These modifications/corrections are listed below and discussed in detail in the sections to follow: (1) derivation of stronger necessary conditions for identifying potential IVs; (2) modification of the procedure for drawing ESDG branches and identification of the conditions under which nodal balance on ESDG cannot be applied; (3) detailed analysis of the cases in which the causal redundant equation derived by nodal balance on the ESDG would potentially help in eliminating spurious solutions; (4) demonstration that CR is implicitly handled by original AEs through appropriate changes in perfect matching while converting a DE to the corresponding AE. 5.2. Stronger Conditions for IR. A system variable xi would show IR with respect to a local exogenous variable el if the signed principal minor corresponding to a directed acyclic path from el to xi is negative. The signed determinant of a matrix H, SD(H), can be written in terms of cycles passing through a variable xi (not to be confused with the variable xi mentioned on the first line above) as follows.29

(-1)r|H| ) -hii(-1)r-1|Qi| -

hjihij(-1)r-2|Qij| ∑ j*i

hjihikhkj(-1)r-3|Qijk| ... ∑ ∑ k*i,j j*i ...∑hjihik...hpj(-1)-1|Qijk,...,p| ∑ p* j*i ...∑hjihik...hpqhqj|Qijk,...,pq| ∑ q* j*i

stronger necessary conditions for xk to show IR with respect to the local exogenous variable el. (1) xk should be located in a negative feedback cycle. (2) The signed determinant of the complementary subsystem to at least one of the directed paths from el to xk should be negative. The recursive necessary condition for condition 2 is that the complementary subsystem contains at least one (i) positive cycle C+ with APSM(C+) > 0 or (ii) negative cycle C- with APSM(C-) < 0. The recursive necessary condition for positive APSM(C) in condition i is that MC should contain at least one (a) positive cycle C+ with APSM(C+) < 0 or (b) negative cycle C- with APSM(C-) > 0. The conditions i, ii, a, and b are recursively used on further subsystems until it terminates (the size of the subsystem becomes 1 or 0). Proof. In eq 48, at least one term (excluding the negative sign before them) on the right-hand side should be positive. This can be achieved in two ways: (i*) the cycle is positive and the APSM is also positive, or (ii*) the cycle is negative and the APSM is also negative. In condition i*, the weaker condition is satisfied at this level itself. In condition ii*, we again use a similar argument for APSM(C). If condition i* is satisfied at any level in which APSM(C+) has order 3 × 3 or higher, then the weaker condition is satisfied. Otherwise, at the next level of decomposition [expansion of APSM(C-)], the order of submatrices would be 2 × 2 or less and the weaker conditions for matrices of order 1 and 2 are easily satisfied.29 Pictorially, in the worst case, the decomposition of the negative APSM looks like the following:

(48)

where Qij is a principal minor of H obtained by removing ith and jth rows and columns from H, etc. By definition, |Qijk,...,pq| is 1 because it has dimension 0 × 0. Also note that this expansion is true for any node i. The necessary condition for (-1)r|H| to be negative is that one of the cycles (including self-cycles) in H should have a positive sign. However, this is not very obvious from the expansion in eq 48. Also, we refer to these conditions as weaker necessary conditions because we have derived stronger conditions later in this section. Oyeleye and Kramer29 have proved the weaker conditions using mathematical induction, but we show a different proof here to motivate the derivation of stronger conditions. The proof of weaker and stronger conditions follows the definition and claim presented below. Definition 5. Associated Principal Signed Minor (APSM). Let A be a square matrix and X be the set of all nodes. Let XC be the set of nodes on a cycle C in A. The APSM of C is the signed determinant of the square matrix MC (formed by eliminating all of the rows and columns from A corresponding to XC).

APSM(C) ) (-1)r|MC| where r is the size of MC and |MC| is the determinant of MC. Claim 6. The necessary conditions presented by Oyeleye and Kramer14 are weak. Given below are the

where C- refers to a negative cycle and Q is a principal minor of order 2 or less. If the positive cycle identified at any of the lower levels passes through a particular system variable xj, then one can write (-1)r|H| in terms of cycles through xj and a positive cycle would be detected at the top level decomposition itself. This completes the proof of weaker conditions. The proof for the stronger conditions continues. For SD(H) to be negative, there must be at least one term on the right-hand side of eq 48 that satisfies condition i* or ii*. If condition ii* is satisfied, then the APSM is required to satisfy condition i* or ii* again. If condition i* has to be satisfied, i.e., a positive cycle is detected, then APSM should be positive and it should satisfy the condition a* or b* given below. The condition for APSM > 0 is as follows: At least one term should be such that (a*) the cycle is positive and the APSM is negative or (b*) the cycle is negative and the APSM is positive. The above procedure is applied recursively until the size of the principal minor is 2 or less, at which point the conditions can be checked by inspection. Note that there is no weaker necessary condition for positive signed determinant. One can check the weaker

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4805

Figure 9. Modified procedure for ESDG construction: modification 2.

necessary conditions for negative SD(H) and, if they are satisfied, then stronger conditions can be checked recursively. Though the procedure for verifying stronger conditions is computationally more expensive as compared to the weaker one, it helps in reducing the set of potential IVs. Also it can be carried out offline while building the qualitative equations for individual units. Next we show an example in which the system matrix (A) contains positive cycles but it never exhibits IR. The weaker necessary conditions predict potential IR, but the stronger conditions do not. The illustrative example is given below. Example 8. Consider the DE system with system matrices

[

]

[]

-7 0 -3 1 A ) -1 1 1 and B ) 0 8.1 -4 1 0

A is a stable matrix. By drawing the SDG for the DE system (not shown here), one can see that weaker conditions predict that all of the three variables, x1, x2, and x3, would potentially show IR with respect to e1. In this case, one would not even know how to draw the ESDG branch. On the other hand, stronger conditions show that only x2 (located in the negative cycle [x1, x2, x3]) and x3 (located in the negative cycle [x1, x3]) can potentially show IR. The stronger necessary conditions are not satisfied by x1. The explanation is simple: one can choose either the self-cycles on x2 and x3 or the cycle between x2 and x3. For the first two positive self-cycles, APSM is not positive, and for the negative cycle, APSM is not negative. In certain cases, when dealing with small systems, it is possible that the stronger necessary conditions might turn out to be the sufficiency conditions as well. The above system is indeed such a system. 5.3. Modified Method for ESDG Construction and Limitations of Nodal Balance on ESDG. Recall that the purpose of drawing ESDG arcs is to provide the negative feedback effect through an acyclic path. Through an example, we will show that the procedure outlined by Oyeleye and Kramer14 for drawing ESDG branches corresponding to IVs does not serve this purpose and hence does not reduce the number of spurious solutions considerably (though it does eliminate some of the spurious solutions).

[

] []

Example 9. Consider the DE system with matrices

-2 1 0 A) -30 0 0

0 0.1 1.6 0 0 0

0 0 -1 -2.5 0 0

1 0 0 -4 3 0

0 0 0 0 1 -1.5

-2 0 0 and B ) 0 0 -2.9

1 0 0 0 0 0

Using the stronger necessary conditions for predicting IR, one can see that x2, x3, and x4 (located in the negative cycle [x1, x2, x3, x4, x5, x6]) show IR with respect to e1. The SDG and ESDG for the above system are shown in parts a and b of Figure 9, respectively. Nodal balance on the ESDG (excluding self-arcs) yields two new (causal) qualitative equations.

[x2] ) [x1]

(49)

[x5] ) [x4] - [x1]

(50)

It turns out that eq 49 is not correct. One can easily see that eq 49 in combination with the AE corresponding to row 2 of A and B predicts X ) 0 as the steady-state solution for nonzero e1, which is wrong. So, we propose the following modifications in the procedure for applying nodal balance on the nodes of the ESDG. Modification 1. If a node has a positive/zero self-cycle on it, then nodal balance should be applied to such a node only if it also receives an IR/CR arc. When eq 49 is ignored, two solutions (the stable steady-state solution and a spurious solution) are obtained. To allow the feedback effect through the ESDG branch corresponding to IV(s), the required modification is as follows. Modification 2. The ESDG branch should start from the local exogenous variable with respect to which certain system variables are exhibiting IR. With the above modification, the ESDG shown in Figure 9c is obtained and the corresponding causal redundant equation (nodal balance on x5) is shown below.

[x5] ) [x4] - [e1]

(51)

Using the above equation with the original AEs, for e1 ) +, the unique stable solution, viz., X ) [+, -, -, +, -, +]T, is obtained. In the above analysis, x5 and x6 (located on the cycle [x1, x4, x5, x6]) show IR with respect

4806 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 Table 6. Usefulness of Causal Redundant Equations: Comparison aii ) 0

a b a b a b

aii < 0 aii > 0

case i (eq 54)

case ii (eq 55)

[xi] ) 0 f useful no new restriction f not useful eqs 54 and 53 are the same f not useful xi can be eliminated from left-hand side f useful

[xi] ) term 3 + term 4 f useful no new restriction f not useful 0 ) term 3 + term 4 f useful no new restriction f not useful 0 ) term 3 + term 4 f useful useful for xi * 0

cases when a causal redundant equation would be useful. Consider the DE

dxi dt

Figure 10. SDG for the DE system: x2 cannot receive an ESDG arc from e1 (modification 3).

to e1 but there is no downstream variable on the cycle so no IR arc is drawn. Next we present another scenario in which nodal balance on some of the nodes of the ESDG is wrong and hence further modifications are proposed. Example 10. Consider a DE system with system matrices (qualitative form)

A)

[ ]

[]

j*i,xj∈Xi

aijxj +



bikek

(52)

k,ek∈Ei

where Xi and Ei are the sets of system variables and exogenous variables, respectively, in the expression for dxi/dt (see eq 47). The corresponding steady-state qualitative equation is

0 ) [aii][xi] +



j*i,xj∈Xi

[aij][xj] +



[bik][ek] (53)

k,ek∈Ei

Next, ESDG is constructed (modifications 1-3 are used), and let us assume that xi receives IR arcs from the local exogenous variables xj′’s and ek′’s. If xj′ ∈ Xi and ek′ ∈Ei, then nodal balance (if applicable, see modification 3) on xi yields (i) eq 54; otherwise, (ii) eq 55 is obtained.

- + and B ) + + -

The SDG for the above system is shown in Figure 10. One can see that x1 shows IR with respect to e1 so the ESDG branch e1 f x2 would have a positive sign. However, the original SDG already contains a negative arc e1 f x2. At this point, one does not know whether to consider the original arc or the ESDG branch. If the original arc is considered, then nodal balance cannot be applied on x2 (modification 1). If the ESDG branch is considered, the resulting system of qualitative equations (including the causal redundant equation obtained by nodal balance on x2) generates the incomplete solution X ) [?, +]T. Numerical simulation for stable A with the above structure generates the solution X ) [?, +]T ∪ [+, ?]T. This set of solutions is also obtained when no ESDG branch is drawn to x2 (and hence no nodal balance is applied). Hence, we propose the following modification. Modification 3. If, on the basis of IR, an ESDG arc is expected from node el to node xj with an appropriate sign and the original SDG contains an arc from the node el to the node xj with the opposite sign, then the ESDG branch is not drawn to the node xj and, hence, nodal balance is not applied on the node xj (modification 1). 5.4. Analysis of the Usefulness of Causal Redundant Equations. The causal redundant equations derived by applying nodal balance on ESDG (with modifications proposed in the previous section) are not always useful in eliminating spurious solutions. For example, if a node with a negative self-cycle receives an ESDG branch due to IR, then the resulting causal redundant equation is merely an expansion of the corresponding original qualitative equation and hence does not impose any additional restriction. In this section, detailed analysis is performed to identify the



) aiixi +

[xi] ) [xi] )



j*i,xj∈Xi

[aij][xj] +



k,ek∈Ei

[bik][ek]



[bik][ek] +

xj′, xi) [xj′] +

sign(IR: ∑ e

[aij][xj] +

sign(IR: ∑ x j′



j*i,xj∈Xi

k,ek∈Ei

(54)

ek′, xi) [ek′] (55)

k′

where sign(IR: xj′, xi) is the sign of the directed arc (xj′ f xi) due to IR. In both the cases i and ii, there are three possibilities: aii ) 0, aii < 0, or aii > 0. Further, eq 53 can be satisfied in two ways: (a) each term provides zero sign individually or (b) there is at least one term providing + contribution to the equation and at least one term providing - contribution to the equation. Thus, there are 12 different scenarios. Table 6 discusses situations in which the causal redundant equations might be potentially useful. In this table, “no new restriction” means that eq 54 or eq 55 contains some additional term(s) as compared to eq 53 without changing the contribution of various terms present in eq 53. In most of the cases the node receiving an IR arc has a positive self-cycle on it, and in such cases as shown in the last row of Table 6, redundant equations are potentially useful. 5.5. Analysis of the CR. In this section, our focus is on analyzing systems exhibiting CR because of the presence of zero self-cycle(s) (that allows one to derive sufficiency conditions). A good understanding of the elaborate work presented by Oyeleye and Kramer14,29 would be helpful in understanding the analysis presented in this section. It should be noted that the conditions derived by Oyeleye and Kramer14,29 are sufficiency conditions rather than necessary conditions. Necessary conditions merely state that either all of the terms on right-hand side in eq 48 are zero or there is at least one positive term and one negative term. The

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4807

Figure 11. SDGs for showing CR: example 11.

Figure 12. Digraphs as qualitative models.

authors (Oyeleye and Kramer) refer to them as necessary conditions as well only in the context of all of the terms being zero, which makes it tangible in the qualitative domain. Thus, the sufficiency conditions for CR are qualitative in nature, unlike the sufficiency conditions for IR that are quantitative in nature. Perfect matchings are considerably different in a DE system (showing CR) and the corresponding AE system because of the presence of zero self-cycle on at least one node. Let us analyze the additional restrictions imposed by a causal redundant equation corresponding to CR. The analysis is similar to the analysis presented in section 5.4. Let xi receive CR arcs from local exogenous variables xj′’s, and let xi* ∈Xi be the CV node just upstream to xi on the concerned negative cycle. The sufficiency conditions show that a node can receive a CR arc in the ESDG only if it has zero self-cycle (i.e., aii ) 0) on it because, if not, then the sufficiency conditions are satisfied for the APSM of the directed acyclic path from xj′ to xi as well. Thus, the first term in the original qualitative equation (53) becomes zero, and the causal redundant equation is

[xi] )



j*i,xj∈Xi

[aij][xj] +



k,ek∈Ei

[bik][ek] +

sign(CR: ∑ x

xj′, xi) δ([xi*])[xj′] (56)

j′

Similar to the arguments used in section 5.4, one can easily see that eq 56 would impose further restrictions on the solution space only if eq 53 is satisfied by [xj*i] ) [ek] ) 0, and then eq 56 becomes

[xi] )

sign(CR: ∑ x

xj′, xi) [xj′]

(57)

j′

Thus, eq 57 indicates that the causal redundant equations corresponding to CR should be useful in some cases, but without rigorous proof, analyzing the following examples shows the contrary. Example 11. Consider the DE system with the system matrices A and B as follows:

[

-3 0 1 0 A) 1 0 0 -1.25 -1

]

[ ]

1 0 0 and B ) 0 1 0 0 0 1

(58)

For the above stable system, the ESDG for the DE system and the SDG for the corresponding AE system are shown in Figure 11. It can be verified that x1 shows CR with respect to e1 and x1 and x3 show CR with respect to e3. So, the CR arcs are as shown in the ESDG. The SDG for the corresponding AE system shows a consider-

able change in the perfect matching. One can write down the causal redundant equation obtained by applying nodal balance on node x2 in the ESDG and simplify it for [x1] ) 0 (corresponding to nonzero e1 or e3). One can also write the original qualitative equations (all three) and can easily verify that for [x1] ) 0, they also reduce to (just by simple substitution and simplification) the causal redundant equation. Thus, the causal redundant equations corresponding to nodal balance on nodes receiving CR arcs are of no additional use for eliminating spurious solutions. The following modification is proposed. Modification 4. CR is automatically handled by steadystate qualitative equations. Another simple example to show the same fact can be constructed by adding a positive arc from x3 to x2 (i.e., [a23] ) +). In this case, x1 and x3 show CR with respect to e3. One can go through a similar analysis to see that the causal redundant equation is useless. Further, we have not been able to prove or refute that a node could not receive both CR and IR arcs from local exogenous variables. Hence, we leave this situation as a possible exception: If a node receives a CR arc, then nodal balance should be applied to this node only if it also receives IR arc(s). 5.6. Discussion. For systems described by DE, IR/ CR has been analyzed in detail in the sections above. In this section we discuss three issues: (i) IR/CR in DAE systems; (ii) IR with respect to the initial response as defined in section 3.1.1; (iii) does our analysis support the validity of nodal balance on ESDG as proposed by Oyeleye and Kramer?29 (i) IR/CR in DAE Systems. Recall that IR for DE systems is predicted based upon stability conditions and that the expression for the ultimate response of a system in terms of directed acyclic paths and APSM requires explicit qualitative information about the system matrix A. For DAE systems, A is not known explicitly. Hence, IR cannot be predicted for DAE systems in general. The only option is to convert the DAE system into a DE system by solving for the algebraic variables and then substituting for them in the DEs (eliminating the algebraic variables). This approach is not feasible for large systems. Notice that one should not calculate time derivatives of the algebraic variables by differentiating AEs because then the resulting system would be singular (see section 3.3.2) and the stability conditions cannot be applied.29 (ii) IR in View of Our Definition of the Initial Response. The ultimate response of a system variable is the total effect of all of the directed paths multiplied with the APSM. In the previous sections, while predicting IR, we effectively tested the necessary conditions for negative APSM. It is important to realize that

4808 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Figure 13. Example A-2: (a) bipartite graph and (b) perfect matchings.

IR (the presence of positive cycles) discussed so far is different from the term IR as used in the process control literature (which is basically realization of the quantitative sufficiency condition of IR). In this section, we refer to IR as discussed in process control terminology. Let us assume that the initial response of a variable can be predicted without ambiguity (i.e., all of the shortest directed paths have the same sign or arc weights can be used to fix the sign). Then a system variable would potentially exhibit IR if (i) at least one longer directed path has the opposite sign and/or (ii) the system variables exhibit IR as defined by Oyeleye and Kramer.29 The conditions mentioned above are the necessary conditions. Thus, if the above conditions are not satisfied for an exogenous variable and a system variable pair, then one can conclude that the chosen system variable cannot show IR with respect to the exogenous variable under consideration. This result can be helpful in applications such as sensor location for incipient as well as steady-state-based fault diagnosis. (iii) Validity of Nodal Balance on ESDG. The following observations are made: (a) Oyeleye and Kramer14,29 have not proved the validity of nodal balance on ESDG. So, one can question its correctness. (b) The methods for identifying IR/CR have been developed to deal with scenarios in which common intuition-based analysis fails; viz., simple propagation is valid for systems containing only negative cycles (or self-cycles). Hence, (without proof) anything else that is based upon intuition also remains doubtful. Further, (c) in the previous section we have shown that the methodology for ESDG construction presented by them requires several modifications. 6. Conclusions In this paper, algorithms for the systematic development of digraphs for DE, AE, and DAE systems have been discussed. Methodologies for digraph analysis have been presented, and the correctness of various algorithms with respect to the corresponding analysis methodology has been established. It is shown that the initial response for DE systems can be predicted by propagating through the shortest paths in the corresponding SDG. Equivalence between steady-state qualitative equations and SDG for AE systems has been established. It is also shown that propagation through the SDG for AE systems is correct if it has only one perfect matching. Analysis of DAE systems is difficult, in general, though great simplification is achieved if the DAE system has only one perfect matching. A detailed discussion has been presented on the use of noncausal redundant equations and causal redundant equations to eliminate spurious solutions in steady-state analysis. IR and CR have been discussed in great detail. A number of examples are provided to clarify and explain the various ideas introduced in this paper.

Acknowledgment Part of this work was carried out at Indian Institute of Technology Bombay, Mumbai, India. Appendix. Definitions 1. Digraphs as Qualitative Models (Figure 12). Digraphs provide a very efficient way of representing qualitative models graphically. A digraph consists of nodes representing variables and directed arcs between nodes representing the interaction among variables. An arc from a node A to another node B means that the variable represented by node A (the initial node) affects the variable represented by node B (the terminal node). There are mainly three kinds of nodes in a digraph representing a chemical process. They are discussed below. Exogenous Variables or Disturbance Variables. These are the variables (denoted by ∀ ei ∈ E) that are not affected by any other variable. They represent disturbance and fault (parametric as well as structural) variables and can change independently. Thus, there are no arcs incident on them. System Variables. These variables (denoted by ∀ xj ∈ X) get affected by exogenous variables and affect each other. So, they have both input and output arcs associated with them. These are often called state variables also. Output Variables. These are the variables (denoted by ∀ yk ∈ Y) that do not affect any other variable, and thus there are no output arcs from these variables. In this paper, output variables have been represented through system variables. An illustrative example follows. Example A-1. Consider flow of a gas through a pipe. The temperature of the gas inside the pipe (T) is assumed to be the same as that of the air outside the pipe (T0). Viscosity µ and density F of the gas depend on its temperature. The cross-sectional area of the pipe, A, is constant. The flow rate of the gas affects the velocity of the gas, which in turn affects the Reynolds number and hence the pressure drop over a certain length (say L; assumed to be constant) of the pipe. This behavior can be represented by equations given below and through the digraph shown in Figure 12. Exogenous variables are T0 and Q. System variables are T, F, µ, and v. ∆P is an output variable.

T ) T0

(A-1)

F ) F(T)

(A-2)

µ ) µ(T)

(A-3)

v ) v(Q,A)

(A-4)

∆P ) ∆P(v,F,µ,L)

(A-5)

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4809

This example consisted of AEs that are encountered during steady-state analysis. 2. Bipartite Graphs and Perfect Matching (Figure 13). Example A-2. Consider the following equations:

f1(e1,x1,x2,x3) ) 0

(A-6)

f2(e2,x2,x3) ) 0

(A-7)

f3(e1,e2,x1,x2) ) 0

(A-8)

Bipartite Graph. A bipartite graph is a graph in which the set of nodes can be divided into two disjoint subsets such that arcs are between nodes of the different subsets only. Let one subset consist of nodes representing the equations and the other subset consist of nodes representing the system variables. We draw an arc from the eq A-6 node to x1, x2, and x3 because eq A-6 contains these variables, and so on for each equation. Thus, we get the bipartite graph shown in Figure 13a. When we delete some arcs such that each equation is matched with one system variable and no equation or system variable is left unmatched, then we get a perfect matching. There can be more than one perfect matching. For the above set of equations, three perfect matchings are obtained, as shown in Figure 13b. Directed Path. A directed path from node A to node B in a digraph is an alternating sequence of nodes and directed arcs of the digraph such that the first and last nodes in the sequence are nodes A and B, respectively. SCC, Maximal SCC (MSCC), and DAG. A subset of a digraphs (consisting of a subset of nodes and arcs of the digraph) is called SCC if every node of this subset can be reached from every other node of this subset. Thus, there is at least one directed path between any two nodes of this subset. A MSCC in a digraph is a node or SCC with no input arcs to it. MSCC is sometimes also called root node. Each SCC can be replaced by a single node called supernode to get DAG. There are no directed cycles in a DAG because if so happens, then there will be again a SCC that can be further replaced by a single supernode. Signed Digraph (SDG). When a sign is put on the nodes and arcs of a DG, it is called a SDG. The sign of a node represents its qualitative state. The values 0, +, and - are assigned to the nodes with normal, high, and low values, respectively. The sign of an arc shows the direction of the effect. A + (-) sign on the arc from node A to node B means that node B would change in the same (opposite) direction as node A. Nomenclature A, A, A′ ) cross-sectional area of the pipe, reactant A, and matrix B, B ) reactant or product B, outlet flow rate from effects 1 and 2, matrix C ) cycle, matrix E ) set of exogenous variables F, F ) volumetric flow rate of reactant, vector-valued function H ) matrix I ) identity matrix L ) length of the pipe P ) product, pressure Q, Q ) flow rate, principle minor S1 ) steam flow rate into effect 1 of a double-effect evaporator T ) temperature

V ) reactor volume W ) holdup X ) set of system variables Y ) set of output variables a2 ) Arrhenius parameter c ) concentration cP ) specific heat e ) exogenous variable f ) molar feed flow rate h ) heat-transfer coefficient, enthalpy, elements of matrix H kr, k2 ) rate constants n ) number of system variables, number of variables in an equation q ) volumetric flow rate r ) size of a square matrix s ) surface area for heat transfer sign(. f .) ) sign of the arc from the first node to the second node t ) time v ) velocity, molar flow rate of the vapor w ) molar flow rate of the bottom product x ) system variable, mole fraction [.] ) qualitative value of the argument First Subscripts A ) reactant B ) reactant or product F ) feed P ) product d ) top product f ) feed i ) inlet j ) jacket i, j, j′, k, l ) indices for variables, equations, and subsets m ) index for variables w ) bottom product 0, 0 ) surroundings, steady state 1, 2 ) subset, effects 1 and 2 Second Subscripts 1, 2 ) more volatile and less volatile components Superscripts -, + ) sign of a cycle Greek Symbols ∆, ∆H ) difference, heat of reaction δ ) Dirac delta function µ ) viscosity φ ) null set F ) density

Literature Cited (1) Rich, S. H.; Venkatasubramanian, V. Model-based reasoning in diagnostic expert systems for chemical process plants. Comput. Chem. Eng. 1987, 11, 111-122. (2) Finch, F. Automated fault diagnosis of chemical process plants using model-based reasoning. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 1989. (3) Wilcox, N.; Himmelblau, D. The possible cause-effect graphs (PSEG) model for fault diagnosissI. Methodology. Comput. Chem. Eng. 1994, 18 (2), 103-116. (4) Wilcox, N.; Himmelblau, D. The possible cause-effect graphs (PSEG) model for fault diagnosissII. Applications. Comput. Chem. Eng. 1994, 18 (2), 117-127. (5) Venkatasubramanian, V.; Vaidhyanathan, R. A knowledge based framework for automating hazop analysis. AIChE J. 1994, 40 (3), 497-505.

4810 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 (6) Venkatasubramanian, V.; Rengaswamy, R.; Kavuri, S. A review of process fault detection and diagnosis. Part II: Qualitative models and search strategies. Comput. Chem. Eng. 2003, 27 (3), 313-326. (7) Iri, M.; Aoki, K.; O’Shima, E.; Matsuyama, H. An algorithm for diagnosis of system failures in the chemical process. Comput. Chem. Eng. 1979, 3 (1-4), 489-493. (8) Umeda, T.; Kuriyama, T.; O’Shima, E.; Matsuyama, H. A graphical approach to cause and effect analysis of chemical processing systems. Chem. Eng. Sci. 1980, 35 (12), 2379-2388. (9) Shiozaki, J.; Matsuyama, H.; O’Shima, E.; Iri, M. An improved algorithm for diagnosis of system failures in the chemical process. Comput. Chem. Eng. 1985, 9 (3), 285-293. (10) Tsuge, Y.; Shiozaki, J.; Matsuyama, H.; O’Shima, E. Fault diagnosis algorithms based on signed directed graphs and its modifications. Ind. Chem. Eng., Symp. Ser. 1985, 92, 133. (11) Kokawa, M.; Shigai, S. Failure propagating simulation and nonfailure paths search in network systems. Automatica 1982, 18, 335. (12) Kokawa, M.; Satoshi, M.; Shigai, S. Fault location using digraph and inverse direction search with application. Automatica 1983, 19 (6), 729-735. (13) Kramer, M.; Palowitch, B. Rule-based approach to fault diagnosis using the signed directed graph. AIChE J. 1987, 33 (7), 1067-1078. (14) Oyeleye, O.; Kramer, M. Qualitative simulation of chemical process systems: Steady-state analysis. AIChE J. 1988, 34 (9), 1441-1454. (15) Chang, C.; Yu, C. On-line fault diagnosis using the signed directed graph. Ind. Eng. Chem. Res. 1990, 29 (7), 1290-1299. (16) Vaidhyanathan, R.; Venkatasubramanian, V. Digraphbased models for automated hazop analysis. Reliab. Eng. Syst. Safe. 1995, 50 (1), 33-49. (17) Mylaraswamy, D.; Kavuri, S.; Venkatasubramanian, V. A framework for automated development of causal models for fault diagnosis. AIChE Annual Meeting 232g, Miami, 1994. (18) Bhushan, M.; Rengaswamy, R. Design of sensor location based on various fault diagnostic observability and reliability criteria. Comput. Chem. Eng. 2000, 24, 735-741. (19) Raghuraj, R.; Bhushan, M.; Rengaswamy, R. Location of sensors in complex chemical plants based on fault diagnostic observability criteria. AIChE J. 1999, 45, 310-322. (20) Bhushan, M.; Rengaswamy, R. Design of sensor network based on the signed directed graph of the process for efficient fault diagnosis. Ind. Eng. Chem. Res. 2000b, 39, 999-1019. (21) Bhushan, M.; Rengaswamy, R. Comprehensive design of a sensor network for chemical plants based on various diagnosability and reliability criteria: 1. Framework. Ind. Eng. Chem. Res. 2002, 41 (7), 1826-1839. (22) Bhushan, M.; Rengaswamy, R. Comprehensive design of a sensor network for chemical plants based on various diagnosability and reliability criteria: 2. Applications. Ind. Eng. Chem. Res. 2002, 41 (7), 1840-1860. (23) Palmer, C.; Chung, P. Creating signed directed graph models for process plants. Ind. Eng. Chem. Res. 2000, 39, 25482558. (24) Chen, J.; Howell, J. A self-validating control system based approach to plant fault detection and diagnosis. Comput. Chem. Eng. 2001, 25, 337-358. (25) Mah, R. Chemical Process Structures and Information Flows; Butterworth Publishers: Boston, 1990. (26) Oyeleye, O.; Finch, F.; Kramer, M. Qualitative modeling and fault diagnosis of processes by MIDAS; Technical Report

LIPSE-89-061; Department of Chemical Engineering, MIT: Cambridge, MA, 1989. (27) Newell, R.; Fisher, D. Model development, reduction, and experimental evaluation for an evaporator. Ind. Eng. Chem. Process Des. Dev. 1972, 11 (2), 213-221. (28) Hall, P. On representatives of subsets. J. London Math. Soc. 1935, 10, 26-30. (29) Oyeleye, O.; Kramer, M. Supplementary Material to qualitative simulation of chemical process systems: Steady-state analysis; Technical Report LIPSE-88-036; Department of Chemical Engineering, MIT: Cambridge, MA, 1988. (30) Basset, L.; Maybee, J.; Quirk, J. Qualitative economics and the scope of the correspondence principle. Econometrica 1968, 36 (3/4), 544-563. (31) Kuipers, B. Qualitative simulation. Artif. Intell. 1986, 29 (3), 289-338. (32) Williams, B. C. Doing time: putting qualitative reasoning on firmer ground. Proceedings of the 5th National Conference on Artificial Intelligence; MIT Press: Cambridge, MA, AAAI 86; Philadelphia, PA, Aug 11-15, 1986; Vol. 1, pp 105-112. (33) Sacks, E. Qualitative analysis of piecewise linear approximation. Int. J. Artif. Intell. Eng. 1988, 3 (3), 151-155. (34) Struss, P. Global filters for qualitative behaviors. Proceedings of the 7th National Conference on Artificial Intelligence; MIT Press: Cambridge, MA, AAAI 88; Saint Paul, MN, Aug 21-26, 1988; pp 275-280. (35) Kuipers, B.; Chiu, C.; Molle, D. T.; Throop, D. Higher-order derivative constraints in qualitative simulation. Artif. Intell. 1991, 51, 343-379. (36) Shen, Q.; Leitch, R. R. Combining qualitative simulation and fuzzy sets. In Recent advances in qualitative physics; Faltings, B., Struss, P., Eds.; MIT Press: Cambridge, MA, 1992; pp 83100. (37) Faltings, B., Struss, P., Eds. Recent advances in qualitative physics; MIT Press: Cambridge, MA, 1992. (38) Bonarini, A.; Bontempi, G. A qualitative simulation approach for fuzzy dynamic models. ACM Trans. Modeling Comput. Simul. 1994, 4 (4), 285-313. (39) Heller, U.; Struss, P. Transformation of qualitative dynamic modelssApplication in hydroecology. 10th International Workshop on Qualitative Reasoning (QR-96), MIT Press: Cambridge, MA, Fallen Leaf Lake, CA, May 21-24, 1996; pp 83-92. (40) Berleant, D.; Kuipers, B. Qualitative and quantitative simulation: Bridging the gap. Artif. Intell. 1997, 95 (2), 215-255. (41) Clancy, D. J. Solving complexity and ambiguity problems in qualitative simulation. Ph.D. Thesis, University of Texas at Austin, Austin, TX, 1997. (42) Keller, U.; Wyatt, T. K.; Leitch, R. R. FRenSi: A fuzzy qualitative simulation method. Workshop on Applications of Interval Analysis to Systems and Control (with special emphasis on recent advances in Modal Interval Analysis), Vehı´, J., Sainz, M. A., Eds.; University of Girona Press: Girona, Spain, 1999; pp 305-313. (43) Kay, H.; Rinner, B.; Kuipers, B. J. Semi-quantitative system identification. Artif. Intell. 2000, 119 (1-2), 103-140.

Received for review August 20, 2002 Revised manuscript received June 10, 2003 Accepted July 7, 2003 IE020644A