3.1.1 The learning strategy of GOLEM

For learning models with GOLEM, QSIM is first formulated in first order logic form, that is, several Horn Clauses, and then the qualitative constraints in QSIM are converted into ground facts, which are the forms that GOLEM can accept. The learning task becomes to construct a model in terms of Horn Clauses that can cover the positive examples obtained from given system behaviours and exclude the hand-generated negative examples:

(2)

GOLEM employs a bottom-up search strategy by iteratively generalizing a most specific clause. For each iteration, first a certain number of clauses are constructed by relative least general generalization (Plotkin, Reference Plotkin1971), among which the clause that can cover the maximum number of the positive examples and exclude all the negative examples will be chosen to be processed for the next iteration (see details in Muggleton & Feng (Reference Muggleton and Feng1990)). Hence, GOLEM is not a complete algorithm in the sense that it uses this greedy hill-climbing heuristic search, and cannot guarantee to find the target model if one exists.

3.1.2 Evaluation of GOLEM for learning QDE models

The validity of GOLEM was tested by learning one of the most commonly used benchmark problems in the QR community: the U-Tube system, which is described in Appendix AFootnote 2. It was reported that on learning the U-Tube system, given a behaviour composed of four states, GOLEM can successfully construct two QDE models in terms of Horn Clauses, both of which are equivalent to the original U-Tube model. Furthermore, GOLEM is capable of finding hidden variables by borrowing the name of available variables (atoms in the clauses).

However, the learning result is based not only on four positive examples, but also six additional hand-generated ‘negative examples’. As mentioned before, the use of negative data is not reasonable in practice, because it is not possible to obtain them unless all the possible behaviours of the system are exactly known. In addition, GOLEM needs a large number of ground facts as background knowledge to support its learning, which inevitably limits its scalability. For learning the U-Tube, there are still 5408 ground facts that have to be generated even after making the following simplifications: (i) the corresponding valuesFootnote 3 are assumed empty; (ii) the consistency of infinite values is ignored; and (iii) multiplication functions are taken out of the constraint repertoire. For dynamic systems that involve long lists of landmark values, the number of ground facts will increase dramatically.

GOLEM demonstrates that ILP is suitable for QML tasks. The idea of employing ILP in QML is adopted by many other later special-purpose QML systems, which explicitly or implicitly employ ILP technology.

3.2.1 Data pre-processing

In contrast to GOLEM, which assumes that the qualitative states are already available, the updated version of GENMODEL possesses the ability to deal with raw numerical data by integrating a Q2Q component. The Q2Q component is made up of two subcomponents: the Front-end Processing module and the Segmenter. The signals pass through the Front-end Processing module, which is composed of several filters, then they are transformed to smooth signals with better mathematical properties. The derivative signals are also obtained by a differentiator. The Segmenter will further divide the smooth signals and the derivative signals, resulting in a complete set of qualitative states.

3.2.2 Model space generation

The model space generation component of GENMODEL is intuitive: the input is a set of qualitative variables (system functions), their dimension information (units), and landmark values. The generation module will exhaustively generate all the possible combinations of the given qualitative constraints and variables, then all the dimensionally inconsistent constraints will be taken out of the constraint set. Under the closed world assumption, the generated model space is complete, which means that it includes all possible models.

3.2.3 Learning strategy

The learning strategy of GENMODEL is straightforward. For each qualitative state obtained from the Q2Q component, all the constraints that are inconsistent with this state will be deleted. After checking all the states, the constraint set is further examined by the redundancy check, which removes the redundant constraints. The redundant constraints are those describing the same relation; for example, M ⁺(A, B) and M ⁻(B, A) are redundant, and one of them should be removed. Finally, the rest of the constraints construct a most specific model that can cover all the qualitative states. To deal with noise, a fault tolerance mechanism is added; each constraint is assigned a counter, recording the number of qualitative states that this constraint is inconsistent with. Only a constraint whose counter exceeds a predefined number can be deleted.

3.2.4 Experimental results

On learning the U-Tube, the early version of GENMODEL (Coiera, Reference Coiera1989a) produced 14 constraints from two behaviours that include six states. Other four physical systems described by Kuipers & Kassirer (Reference Kuipers and Kassirer1984) and Kuipers (Reference Kuipers1986), including the bathtub and the spring systems, were also tested by the same version of GENMODEL. Experimental results showed that this version of GENMODEL was able to learn these models, but GENMODEL tended to generate overconstrained models when the given data were insufficient. Furthermore, the updated version of GENMODEL (Hau & Coiera, Reference Hau and Coiera1993) was successfully applied to a real-world application, which will be detailed in section 6.2.

3.2.5 Evaluation

One of the advantages of GENMODEL is that it does not need negative examples, and the dimensional information can serve as a form of directed negative example generation. It also has the following features: (i) It has the ability to handle real-valued experimental data; (ii) it can deal with noisy data to some extent because GENMODEL introduces the fault tolerance mechanism described in section 3.2.3; (iii) unlike GOLEM, it does not need ground facts; and (iv) it is a complete algorithm, which means that it can guarantee to find the target model if one exists.

The limitations of GENMODEL are as follows: (i) it cannot introduce hidden variables. If not all the system variables are identified, GENMODEL will terminate at a ‘shallow model’. (ii) As it tends to find the most specific models, the models obtained are usually overconstrained, that is, it contains more constraints than necessary to fully characterize the dynamic system being modelled.

GENMODEL’s primitive idea of generate-and-test is also adopted by later QML systems. These learning systems also added other techniques based on other considerations. For instance, the constraint determination stage of QSI (Say & Kuru, Reference Say and Kuru1996) is basically the same as GENMODEL; MISQ (Richards et al., Reference Richards, Kraan and Kuipers1992) also uses the same method as GENMODEL to generate its initial constraint set.

3.3.1 Data pre-processing

This component is called the Conversion of qualitative data phase in MISQ. This module can take as input the data from a high-resolution sensor, or quantitative behaviours directly generated by hand, which are composed of those time points where some variables reach extreme values or zero. If the inputs are sensor data, they will be converted to quantitative behaviours by extracting the significant data points. The quantitative behaviours will further be transformed to qualitative behaviours by constructing the landmark values and estimating the derivatives. This phase will be bypassed if the input data are qualitative behaviours.

3.3.2 Model space generation

The Model Space Generation (Constraint generation phase in MISQ) is basically the same as GENMODEL, that is, it first generate all the possible combinations of qualitative constraint type and variables, then use dimensional information to filter out the inconsistent constraints.

3.3.3 Learning strategy

MISQ will generate one and only one most specific model if the input qualitative behaviours to this module satisfy the following two conditions: (i) These qualitative behaviours are consistent. A set of behaviours is consistent if these behaviours can be achieved by the same real dynamic system. Note that MISQ does not require complete behaviours. (A set of behaviours is complete if it includes and only includes all possible behaviours that a real dynamic system can achieve.) Hence, this is a more relaxed condition than that in GENMODEL, which requires that the input qualitative behaviours must be complete in order to obtain the right model. (ii) The descriptions of these behaviours are complete, in the sense that there are no hidden variables in the system, and there are also no missing values and dimensions for all the qualitative variables. Considering that condition (i) only requires the consistency of the input qualitative behaviours, the input can be a subset of the complete behaviours, in this case MISQ may produce an overconstrained model.

MISQ can also learn models when the above two conditions are not satisfied, that is, (a) there are hidden variables and (b) there are missing qualitative values or dimensions. In the case of (a), MISQ will use the relational path-finding algorithm (Richards & Mooney, Reference Richards and Mooney1992) to find the hidden variables. In the case of (b), an inconsistent constraint set may be generated. All the subsets of this constraint set construct a QDE model search space. MISQ will search in this space by employing some heuristics, such as ‘the subset must contain at least one derivative constraint’ (the system being modelled must be dynamic) and ‘form a connected graph’ (no disjoint models). These heuristics are also adopted in a later system, ILP-QSI (Coghill et al., Reference Coghill, Garrett and King2004, Reference Coghill, Srinivasan and King2008), in which the well-posed models are formally and completely defined to include more properties that a QDE model should have.

3.3.4 Experimental results

On learning the U-tube, MISQ inferred the right model from the only one behaviour, which was the same used in GOLEM as positive examples. Experiments on learning the cascaded tanks (described in Appendix C) were also performed when given data were at different resolutions (qualitative, quantitative, and high-resolution data), or had hidden variables. Results from these experiments showed that MISQ could produce the righ-cascaded tanks model if the dimensional information was complete, and if there were missing variables, MISQ could also infer the right model by finding the hidden variables through the relation path-finding approach.

3.3.5 Evaluation

MISQ goes further than GENMODEL by possessing the ability to process incomplete information and deal with hidden variables. It can work with both quantitative and qualitative data. One of the characteristics of MISQ is that it allows partially specified variables in a qualitative state, which is useful when measuring such variables is not very easy or some points of measurements for them are missing. All of the above make MISQ an advanced learning system that holds an important position in the field of QML.

The limitations of MISQ are as follows: (i) Like GENMODEL, even when complete behaviours are provided, MISQ may produce overconstrained models because it tends to generate the most specific models. (ii) When dealing with incomplete information (a), the heuristics used to reduce the model search space are not completely and formally defined; so additional heuristics can be added, such as causal ordering (Iwasaki & Simon, Reference Iwasaki and Simon1986), which is also used in a later learning system ILP-QSI (Coghill et al., Reference Coghill, Garrett and King2004, Reference Coghill, Srinivasan and King2008). This can help it overcome the generation of overconstrained models; (b) no further systematic experiments were performed to support its learning ability, although some examples such as learning the cascaded tanks under a few types of input were demonstrated. (iii) It is not apparent from the descriptions or experimental results of MISQ whether or not it can handle noisy data. The Correctness and Uniqueness Theorem presented only applies for complete noise-free data.

3.3.6 MISQ-RT (MISQ for region transitions)

MISQ-RT (Ramachandran et al., Reference Ramachandran, Mooney and Kuipers1994) can be seen as an extended version of MISQ, and it can deal with hybrid dynamic systems (Alur et al., Reference Alur, Courcoubetis, Halbwachs, Henzinger, Ho, Nicolin, Olivero and Sifakis2000). Briefly, a hybrid dynamic system is composed of several continuous dynamic systems and the switching conditions among them. In the context of qualitative modelling, each of these dynamic systems can be modelled by an individual QDE. Each QDE holds under specific conditions, which are termed operating conditions of the QDE, and expressed by the values of the variables involved. The operating conditions specify a ‘region’ within which this QDE holds, and this region is called the operating region (Kuipers, Reference Kuipers1994). If a switching condition is satisfied, a region transition will occur, which means that the hybrid dynamic system will switch from the operating region of one QDE to that of another.

MISQ-RT fulfils the task of learning hybrid dynamic systems by executing the following four consecutive steps: (i) Detecting region transition points by heuristic methods, such as detecting the discontinuous change of the magnitude and the sign of the derivative. Then, the behaviours are divided by these transition points into segments. Each segment is an operating region, or a subset of an operating region. (ii) Learning QDE models within each segment by MISQ. (iii) Identifying the operating conditions for each region. The operating condition for an operating region is initially set as the conjunction of the range of all the qualitative values for each variable within this region. (iv) Unifying the regions to get a compact description of the model. In the first three steps, for each behaviour segment, a corresponding QDE is obtained. However, some of the segments may have the same QDE model. So it is necessary to unify these segments. Two heuristics are used: (a) The identical constraints heuristic will try to unify two regions with the same QDE model, and the operating condition for the newly unified region is updated according to those for the original regions; (b) the identical operating conditions heuristic will try to unify two regions with the same operating conditions. The QDE model in the new region is the intersection of the models in the original regions.

As far as we are aware, MISQ-RT is the only QML system that can deal with hybrid dynamic systems.

3.4.1 Data pre-processing

There is no explicit description of a Q2Q component in QSI. However, QSI introduces a not fully implemented qualitative noise filter. This filter tries to smooth the qualitative behaviours by replacing those qualitative states with unnecessarily noisy landmarks with a more general qualitative state. A parameter is specified to indicate the sensitivity of the filter, that is, the length of the qualitative state sequence to be smoothed. If within this given length, a sequence is the so-called tooth one, that is, for one or more than one step a variable increases (or decreases) first, then becomes steady, and finally decreases (or increases), the filter will replace this sequence with a single qualitative state.

3.4.2 Initial search space generation

The initial search space generation is implemented by a subcomponent named constraint determination, which is also reused by the later search process. Similar to GENMODEL, constraint determination tries to combine all the input variables and constraint repertoire to generate all the possible qualitative constraints. In contrast to GENMODEL, first, as QSI can start from a subset of all the system variables and postulate the missing ones, this subcomponent can take even a small fraction of all the system variables as input. Second, as QSI assumes that the input states are complete and noise free, this subcomponent will remove those generated constraints, which are inconsistent with any of the input states. Third, the newly generated qualitative constraints, which are the consequences of the existing constraints, are also ignored.

3.4.3 Learning strategy

As mentioned above, QSI employs a unique iterative search approach over the model space: it will first try to construct the model using only the given variables. If the given variables are only part of all the variables in the system, the constructed intermediate model will be more general than the target model, which means that it will output more behaviours than expected. This model is termed a ‘shallow’ model in QSI. Then the algorithm will try to find a more specialized model based on this ‘shallow’ model, which can cover a few other behaviours, which are not included in the input qualitative states, as possible. In order to achieve this, the following three steps are executed iteratively until a satisfactory model is found:

Step 1: The current intermediate model undergoes the model depth test. The model depth test employs a slightly modified QSIM simulator to verify the current intermediate model. If this intermediate model is a ‘shallow’ model, the QSIM simulator will output more behaviours than those input to QSI. In this case, a model extension process (Step 3) is necessary.

Step 2: If the current model passes the model depth test, a dimensional consistency operation will be performed: Additional M ⁺ constraints are added to the model to eliminate the dimensional inconsistency. In each of these M ⁺ constraints, one argument is a variable in a dimensionally inconsistent constraint, and another argument is a newly introduced variable (the so-called ‘buffer parameter’), which has the same quantity space structure as the other variable in this M ⁺ constraint but has a different dimension. Lastly, QSI will return the final resulting model.

Step 3: If the current model fails to pass the model depth test, the model extension will be performed:

• • Step 3.1: New hidden variables are postulated. This step can be run in several postulation modes to control the computational expense. In the worst case, full postulation mode, which postulates new variables by generating the defining constraintsFootnote 4 from all combinations of the constraint repertoire (except the monotonically increasing and decreasing functions) and available variables, will be employed if no additional knowledge about the model is known. For example, on learning the U-tube, two new variables may be postulated: P_X and P_Y, and their defining constraints are as follows: DERIV(amount 1, P_X), ADD(amount 1, amount 2, P_Y).

• • Step 3.2: The newly postulated variables will be assigned behaviours according to the following two heuristics: the newly assigned behaviour should have the minimum number of qualitative direction changes, and the postulated variables should be constant whenever possible. Thus, a wider and bigger qualitative state set, described by all the variables (including the newly postulated ones), is obtained.

• • Step 3.3: Extended constraint determination will be performed on this qualitative state set. This will return a new model for the next iteration. Again because of computational expense considerations, different search modes are provided; the full search mode will try all the combinations, while half search mode will only try the combinations involved in at least one old variable.

3.4.4 Learning the U-tube as an example

On learning the U-tube, suppose the given data contained the information about only two variables, the amount of liquid in tanks 1 and 2, amounts 1 and 2, a model consist of only one constraint M ⁻(amounts 1 and 2) was first constructed. This model failed to pass the model depth test (Step 1) as it was too ‘shallow’. Then the model extension (Step 3) was executed. If in Step 3 the half search mode was applied and only the derivatives of existing variables were postulated, the QDE model after the execution of Step 3 is shown in Table 2, where P ₁ and P ₂ are newly postulated variables.

Table 2

The U-Tube QDE after the first iteration of QSI

QDE, Qualitative Differential Equation; QSI, Qualitative System Identification.

In the second iteration, the model shown in Table 2 passed the model depth test in Step 1, and the dimensional consistency (Step 2) was performed subsequently, resulting in a final output model, which is dimensionally consistent.

3.4.5 Evaluation

Possessing most of the desired features of a QML system, QSI can be deemed a state-of-the-art QML system. The unique search algorithm makes QSI a distinctive QML system, and also gives it the ability to deal with hidden variables effectively.

The limitations of QSI are as follows: (i) Unlike MISQ, QSI is unable to handle incomplete input behaviours. QSI requires that the input must be correct and complete, otherwise it may output a wrong model. (ii) It cannot deal with quantitative data. The input must be in qualitative states. (iii) The noisy data processing is not fully implemented. Although reasonable, the qualitative noise filter needs further development. (iv) QSI cannot easily incorporate additional constraints upon the model to guide its search, as ILP-QSI does.

3.5.1 Data pre-processing: Q2Q

The core algorithm of ILP-QSI can learn from only qualitative data, but quantitative to qualitative conversion is also provided. The central difference approach (Shoup, Reference Shoup1979) is used to estimate the first and second derivatives of a quantitative variable. The following simple formulae are used:

(3)

(4)

In the above formulae i = 2,3,…,N − 1, and N is the number of input time points. Then a Blackman filter (Blackman & Tukey, Reference Blackman and Tukey1958) is used to smooth the estimated first and second derivatives. Finally, a quantitative variable x is converted to a qualitative variable q = 〈qmag, qdir〉, where qmag is generated from x and qdir is generated from dx/dt. The qualitative derivative of q and q′, is obtained in a similar way, whose qmag and qdir are generated from dx/dt and d ²x/dt ², respectively.

3.5.2 Learning strategy

As mentioned above, the search algorithm of ILP-QSI is implemented within the ALEPH system, and is a variant of the branch-and-bound algorithm (Papadimitriou & Steiglitz, Reference Papadimitriou and Steiglitz1982). The search strategy is detailed as follows:

• • In ILP-QSI, a QDE model is represented as a definite clause.

• • Constraints upon the model, for example, the well-posed model constraints (described in the next section), are assumed to be encoded in the background knowledge.

• • The algorithm starts from an empty set, while maintaining a prioritized queue Active, which is initially empty and later will store the explored models, which are sorted by the cost function.

• • The cost of function of a model is calculated from

(5)

where P(C|B,E) is the Bayesian posterior probability estimate of clause C, given the background knowledge B and positive examples E. The model with the maximal posterior probability is approximated by maximizing the following function (McCreath, Reference McCreath1999):

(6)

where D_H is a prior probability, g is the generality of a model, p is the number of positive examples. These functions and parameters can be calculated in the manner of C-Progol (see Muggleton (Reference Muggleton1996) for details).

• • Successor Function: The unexplored nodes are expanded by a successor function, which can efficiently enumerate the well-posed models. The successor function is defined to return a set of all the well-posed models when the current node is empty and return empty when the current node is not empty. All the returned nodes are added to the branch set.

• • Branching: (i) A node k from the queue active is popped out and evaluated by the above-defined cost of function. (ii) The cost of the best node so far, best, is updated if necessary. (iii) Node k is expanded by the above-mentioned successor function. All the expanded nodes are added to the branch set.

• • Bounding: Prune the branch set according to the cost best and obtain a bounded set, then add the elements in the bounded set to the queue active.

• • The branch-and-bound operation executes iteratively until the queue active is empty or the number of explored nodes exceeds a pre-specified number.

3.5.3 Well-posed models

The concept of the well-posed model was proposed by analyzing the properties of a class of dynamic systems in depth, and it reveals the internal mechanism of these systems; on the other hand, well-posed model constraints can reduce the search space; the branch-and-bound algorithm will only expand and evaluate those nodes, which satisfy the well-posed model constraints. A well-posed model must satisfy the following syntactic constraints:

• • Size: The size of the model must be specified; this is obtained by counting the number of qualitative constraints in the model.

• • Complete: The model must include all the observed variables.

• • Determinate: The model must contain as many qualitative relations as qualitative variables.

• • Language: The number of instances of any qualitative relation in the model must be below some pre-specified limit. Detailed discussion can be found in Camacho (Reference Camacho2000).

The well-posed model must also satisfy the following semantic constraints:

• • Sufficient: The model must adequately explain the observed data.

• • Redundant: There are no redundant relations in the model.

• • Contradictory: There are no contradictory relations in the model.

• • Dimensional: All the constraints in the model must be dimensionally consistent. The details of dimensional analysis can be seen in Bhaskar & Nigam (Reference Bhaskar and Nigam1990).

In addition, the following constraints are preferred but not mandatory:

• • Single: The model must not contain disjoint models, that is, the target model must describe the qualitative behaviours of a dynamic system within the same operating region.

• • Connected: All the hidden variables should appear in at least two qualitative constraints in the model.

• • Causal: The model must be causally ordered (Iwasaki & Simon, Reference Iwasaki and Simon1986; Wiegand, Reference Wiegand1991).

3.5.4 Experimental design and results

Using the above-mentioned search strategy, together with the well-posed model assumption, a series of experiments are tested.

Learning without noise: Noise-free data were provided to the learning system and the experiments were tested by the following three steps: (i) Obtain the complete envisionment (Kuipers, Reference Kuipers1994) of the system with specified values of exogenous variablesFootnote 5. (ii) Power Set Experiments: all the non-empty subsets of the qualitative states in the complete envisionment were treated as training data. If the number of qualitative states was N, the number of experiments was 2^N − 1. For small N, ILP-QSI exhaustively tested on these training data sets. For large N, a fixed number of experiments was randomly sampled from the power set of all the qualitative states. (iii) The learning reliability was measured by the precision of learning from these training data sets. The average precision under different sizes of training data was obtained and visualized.

Learning in the presence of noise: Non-empty subsets of the states in the complete envisionment were replaced by randomly generated noisy qualitative states. All the possible replacements, with the number of 2^N − 1, were tested and again the average precision over different sizes of noisy data was visualized.

Kernel subsets and solution space analysis: With noise-free data, a kernel subset in ILP-QSI is a subset of the complete qualitative states, and is composed of minimal qualitative states required for successful learning (i.e. the target model can be successfully identified in an experiment). The meaning of this is twofold: first, any superset of a kernel subset will also lead to successful learning; second, any subset of a kernel subset will fail to identify the target model. For example, on learning the U-Tube, the kernel subsets are eight pairs of qualitative states. This means that in all the successful experiments, the training data must be a superset of at least one of these eight pairs.

A further investigation by analyzing the solution space of the U-Tube and coupled tanks (described in Appendix B) reveals that the qualitative states in the kernel subset include at least one critical point and these states cover different branches in the envisionment graph (Coghill et al., Reference Coghill, Srinivasan and King2008).

3.5.5 Evaluation

ILP-QSI is an ILP-based QML system, which retains all the characteristics of the earlier QML systems, such as MISQ and QSI. In addition, ILP-QSI has its own new features: (i) the concept of well-posed models is explicitly and completely defined; (ii) it can learn in the presence of noise; (iii) It can avoid finding overconstrained models; and (iv) it can be used to learn complex biological systems. The metabolic pathway of glycolysis was chosen to test the learning ability of ILP-QSI on large systems, and this will be described in section 6.4. All of the above make ILP-QSI a state-of-the-art QML system. In addition to the above-mentioned new features, a novel contribution of ILP-QSI to the QML research is the kernel subset experiments that were performed.

There are, however, a number of ways in which ILP-QSI could be improved: (i) Like many other QML systems, it has not been tested to deal with the dynamic systems with multiple operating regions. (ii) Another problem is the computational limitation of ILP-QSI for learning complex systems. In identifying the pathway of glycolysis, several months of compute time have to be used. This will also be discussed in section 6.4.

3.6.1 Data pre-processing

As QMN directly learns the QDE models from numerical data, the data pre-processing is quite simple; it mainly involves estimating the time derivatives of the system variables. Given the time series points x_i (i = 1,2,…,N), the derivatives are obtained by the following formula from (Bohte, Reference Bohte1991):

(7)

where h is the time interval between two time points. Given the highest order o, QMN will generate all the o derivatives for a variable using the above formula. Then the time series points and their derivatives are directly used by QMN. It should be pointed out that there is no noise processing in the pre-processing.

3.6.2 Learning strategy

The core algorithm of QMN uses the generate-and-test strategy, and is based on GENMODEL, except that it can introduce new variables in a simple manner. (i) First, all the variables and their derivatives are added to an initial set V. (ii) Introduce new variable: Given the model depth d, new variables will be iteratively introduced by combining the variables and their derivatives in V. Four basic arithmetic operations: addition, subtraction, multiplication, and division are used to combine any two elements in V. The newly introduced variables are added into V for the next iteration. (iii) Generate-and-test: Given two tolerance parameters δ and ε, any possible combination of QSIM constraints and variables will be generated and tested. All the generated constraints which are dimensionally consistent are transformed to their equivalent constraints in the form zero(T), which holds when the term T is considered as a constant zero given the tolerances δ and ε. For example, add(X,Y,Z) is converted to zero(Z−X−Y), and if zero(Z−X−Y) holds, the original constraint add(X,Y,Z) also holds. Then the test of the validity of a constraint is converted to test whether zero(T) holds or not. This is achieved by testing the validity of P(|T| > δ) < ε, where P is the probability function, which can be approximated by calculating the percentage of the measurements in which |T|>δ holds over the total number of measurements.

3.6.3 Evaluation

Instead of learning from qualitative data, QMN can directly handle numerical data, and this is achieved by introducing two tolerance parameters and the corresponding test method. The highest order of the system o and model depth d determine the complexity of the search space. Careful selection of these parameters is very important for learning.

Like GENMODEL, QMN tends to generate most specific constraints; this often leads to an overconstrained model. Experiments on the U-Tube and cascaded tanks system (described in Appendix C) showed that not only redundant constraints, but also some unexpected constraints were generated. It was reported that for the U-Tube with d = 2, there were 362 constraints generated. At this point, a further improvement for QMN is required to reduce the number of constraints generated. Another limitation of QMN is that it needs a large amount of high-resolution quantitative data. For learning the U-Tube, 1000 qualitative data points had to be provided to QMN. This restricts its learning ability when only sparse data are available.

3.7.1 Data pre-processing: interval partition and polynomial approximation

The data pre-processing is a significant part of LYQUID. For each of the system variables, given the observed data samples, together with the time points of measurements, the interval partitioning algorithm will try iteratively to partition the time intervals into subintervals, and approximate the variable within these intervals by individual orthogonal polynomials. All the intervals on which the square error of the current polynomial is greater than a specified parameter Tolerance will be dealt with by the interval partitioning algorithm.

When the interval-partitioning algorithm tries to split a time interval, first it temporarily splits this interval into two subintervals from the midpoint, then uses two new polynomials to fit on these two subintervals separately. If the improvement of the precision upon the original interval is greater than the pre-specified parameter improvement, or the improvement of precision on either of the two subintervals is greater than the pre-specified parameter split improvement, this split is accepted; otherwise, the split is rejected and the original interval will be kept. The precision of the polynomial in an interval is calculated by the square error of this polynomial for approximating the variable within the interval over the length of the interval. The newly split intervals will also be recorded for further processing by the partitioning algorithm. The partitioning algorithm will terminate when no further intervals can be split.

3.7.2 Learning strategy

Like QMN and GENMODEL, LYQUID uses the generate-and-test strategy; all the possible combinations of qualitative constraints and variables are tested. Similar to QMN, testing whether a qualitative constraint holds or not is equivalent to testing whether two polynomials are equal, or strictly speaking, their difference is within a tolerance interval. LYQUID uses the following formula to measure the difference between two polynomials p and q on a time interval (a and b):

(8)

If e is less than a specified parameter iTolerance, polynomials p and q are deemed to be equal. iTolerance is the specified tolerance square error per unit time. As in QMN, the parameter o specifies the highest order of a model. The derivatives of a variable are generated by differentiating the corresponding polynomial that represents this variable. As LYQUID evaluates the validity of a qualitative constraint on many time intervals, and this qualitative constraint may not hold on all the time intervals, the ratio R = s/k is calculated by LYQUID, where s is the number of intervals on which the constraint holds, and k is the total number of intervals. Additional parameters have to be set to determine whether the constraint holds or not, and if the ratio R falls into the range determined by these parameters, the qualitative constraint can be considered to hold on all time intervals.

The steps of the LYQUID core algorithm are as follows: (i) For all the observed variables, add the corresponding polynomials and the dimensional information into the set of variables V. (ii) Hidden variable generation: generate the polynomials for the derivatives, sums, and products of the variables in V. Add these newly generated polynomials into V. (iii) Check all the possible combinations of QSIM constraints and the variables in V, adding the consistent constraints into the constraint set C. (A consistent constraint is defined as the constraint that agrees with the data. More specifically, a consistent constraint in LYQUID is the constraint in which the difference between its corresponding polynomials on both sides is within the tolerance interval.) The algorithm will finally return the constraint set C.

3.7.3 Experiments

LYQUID was tested on the U-Tube, cascaded (see Appendix C) and coupled tanks (see Appendix B), the same physical systems used in ILP-QSI. The algorithm was tested on both clean and noisy data samples. In addition, the experiments on the data, which are obtained from irregular sampling frequencies, and the experiments on the missing data were tested.

3.7.4 Evaluation

LYQUID, as a QML system, which can learn models directly from numerical data, uses polynomials as a tool to fulfil the learning task. Because of the nature of polynomial approximation, LYQUID has the following features: (i) it is suitable for dealing with noisy data; (ii) the variables can be sampled at different frequencies; and (iii) It can also work on incomplete data, or even arbitrary sampled data.

A major limitation of LYQUID is that many parameters have to be specified properly by the users. These parameters are even more than those in QMN. Another limitation is that, like GENMODEL and QMN, LYQUID often generates overconstrained models with unnecessary constraints.

3.8.1 Data pre-processing

On the basis of the above assumptions, a set of rules has been established to achieve the following two goals: (i) distinguishing the time intervals and moments from the input sequence; and (ii) estimating the qualitative states. The latter is carried out by estimating the possible qualitative derivatives according to the input time series of qualitative values. Two state estimation rules have been set up, depending on whether the current quantized time is a time interval or a time point. Table 4 shows part of the State Estimation Rule 1, when the current quantized time is a time point. Table 5 is part of the State Estimation Rule 2, when the current quantized time is a time interval.

Table 4

Part of the State Estimation Rule 1 in Abe’s System

Table 5

Part of the State Estimation Rule 2 in Abe’s System

In Table 4, [Q(t_now)] stands for the qualitative value at the current time point; [Q(t_prev)] and [Q(t_next)] stand for the qualitative values on the previous and succeeding time intervals separately. ∂[Q(t_now)] is the estimated derivative. In Table 5, [Q(t_now)] is the qualitative value within the current time interval, so [Q(t_prev)] and [Q(t_next)] are both qualitative values on time points. ∂[Q(t_{now −s})] is the estimated derivative in the interval immediately after the time point t_prev; ∂[Q(t_{now −e})] is the value in the interval immediately before the time point t_next; ∂[Q(t_{now −m})] is the value in the time interval between t_{now −s} and t_{now −e} and Q(t) is assumed to be strictly monotonic in t_{now −m}.

By applying the State Estimation Rules 1 and 2 on the input sequence data, all the possible qualitative states can be obtained.

3.8.2 Model space generation

All the possible qualitative differential equations, which are different from those used in QSIM, are represented by the following expressions:

(9)

where, ⊂ means ‘subset’. For example, [Q]⊂[0+] means [Q] = [0] or [Q] = [+]. ∑ and ∏ means qualitative addition and multiplication, respectively. order is the highest order of the polynomial, and [C_{k 1,…,kn ,l 1,….,ln}] are the coefficients. All possible qualitative differential equations indicated by equation

are tested for validity by all the estimated qualitative states, which are obtained by the two state estimation rules described in section 3.8.1. Finally, the validated constraints that are consistent with the estimated qualitative states constitute the search space.

3.8.3 Learning strategy

The learning task is to find the so-called Simultaneous Differential Equations, which are the conjunction of the qualitative differential equations in the above search space, to construct a model. Abe’s system uses the concept of Entropy Reduction (Cellier, Reference Cellier1991) as a heuristic function to select the appropriate Simultaneous Differential Equations. A bigger value for entropy reduction means a stronger forecasting power. Entropy reduction is defined as:

(10)

where n is the number of variables and k is the number of qualitative states that satisfy the model. When n is fixed, a smaller k will result in a larger entropy reduction H_r. Then the models that can generate fewer states will be chosen as the solution models.

3.8.4 Experimental results

A randomly generated input sequence of qualitative values with two variables is selected as the test example. It is reported that six models are found. All these models are transformed to QSIM models and simulated by QSIM. They can reproduce the same data as the input.

3.8.5 Evaluation

Abe’s QML system is based on several assumptions and particular conditions, and these assumptions and conditions can reduce the complexity of the learning system and simplify the learning task. On the other hand, because the algorithm restricts itself to the above conditions, it will encounter some difficulties when applied to some real applications that do not obey the above conditions. For example, when the given data are sparse, incomplete, or there are hidden variables, the algorithm will not be able to learn the right model.

3.9.1 An incremental decomposition strategy

QML-IL uses a simple incremental decomposition strategy to decrease the computational complexity of ILP-QSI. The learning task is divided into several stages and the decomposition is achieved by domain-specific knowledge. Each stage is provided with the following parameters: (i) the set of variables involved in this stage; (ii) the number of constraints in this stage.

The background knowledge B and examples E are transformed into a sequence 〈B_m,E_m〉, where m = 1,2,…,n, and n is the total number of stages. Each of the elements in this sequence, 〈B_i,E_i〉, is provided to the corresponding stage i. Where E_i contains the observations for the variables involved from stage 0 to stage i, and B_i is the background knowledge that involves the corresponding variables in E_i.

In each stage, say, stage i, given 〈B_i,E_i〉, and H_{i − 1}, which is the set of partial models that are identified in stage (i − 1), the learner L will output the set of partial models H_i. This process can be denoted as the following formula:

(11)

where i = 0,1,….,n, and n is the total number of stages; f is the cost function as described in equation (5). H ₀ is an empty set, and H_n is the set of final models. ρ is the successor function defined as follows: Given each element k in H_{i − 1}, ρ(k,B_i,E_i) will return all models that obey the well-posed model constraints and the number of qualitative constraints in the current stage.

The implementation within the learner L is exactly the same as ILP-QSI, which uses the branch-and-bound search strategy. In addition, before the set of partial models, H_{i − 1,} is submitted to L in stage i, there is an optional operation that performs a limited generalization on each model in H_{i − 1}. This is implemented by disconnecting some or all components in each model, but respecting the original constraints provided in this model. In ILP, this is called removing variable co-references. After the learner L returns a set of models H′_i in stage i, there is another optional operation: randomly selecting a certain number of models from H′_i, and these randomly selected models will compose H_i, which is used for the input of the next stage.

For example, for learning the coupled tanks system described in Appendix B, the learning task can be divided into two stages by assuming that tanks 1 and 2 are relatively separated (domain-specific knowledge): The first stage concerns only the variables related to tank 2, and tries to identify the models that are consistent with the observations; the second stage concerns all the variables and tries to learn the final model based on the models found in stage 1.

3.9.2 Automatic decomposition

In some real-world problems, the decomposition of a dynamic system may not be made appropriately, and thus the automatic decomposition algorithm is required. The decomposition task is actually to determine the partition sequence (S ₀, S ₁, S ₂,….,S_n) on the set of all the given system variables, namely S, where S_i (i = 1, 2,…,n) is a subset of S, and S = S ₀ ⋃ S ₁ ⋃ S ₂… ⋃ S_n.

A randomized local search algorithm, which is similar to GSAT (Selman et al., Reference Selman, Levesque and Mitchell1992), is used to perform the automatic decomposition. The steps are briefly described as follows:

Suppose the first (i − 1) subsets of variables from S for the first (i − 1) stages have been selected: S ₀, S ₁, …, S_{i − 1}; S ₀ = φ; VarsLeft = S, when i = 0; and VarUsed = S ₀ ⋃ S ₁… ⋃ S_{i − 1}. Then for the selection of S_i:

1. 1. First, randomly select a subset V from the unselected variables VarLeft = S-VarUsed for current decomposition.

2. 2. Determine the cost of this selection c by the incremental learner L. c is the cost of the model which has the least cost in H′_i, and H′_i = L(H_{i − 1}, B, E, VarUsed ⋃ V). Updating the current best cost, bestcos,t by c and current selection, VarsSelected, by V, if c is lower than bestcost.

3. 3. BestLocal = V.

4. 4. Construct all the possible ‘local moves’ for BestLocal. This is implemented by removing one variable from V and adding another variable, which is not included in V.

5. 5. Determine the cost of all these moves and selecting V′, which has the least cost c′.

6. 6. Update bestcost and Varselected with c′ and V′, if c′ is lower than bestcost.

7. 7. BestLocal = V′.

8. 8. 4–7 are executed M times.

9. 9. 1–8 are executed R times.

10. 10. S_i = VarsSelected.

In the above steps, M and R are very important parameters. M is the ‘depth’ of the local moves and R is the number of restarts of the randomized procedure. An initial investigation of the influence of these two parameters is also provided: both the coupled tanks system and other three artificial systems are tested by the automatic decomposition algorithm. The experimental result indicates that (i) generally, the algorithm will perform better if R≥M; (ii) when the number of variables in the system or the number of stages is increased, R and M should also be increased.

3.9.3 Experiments and evaluation

Besides the coupled tanks system, the models of different biological systems were also tested. They are listed as follows: (i) predator–prey model; (ii) human lung model; and (iii) glycolysis pathway (the same model is used in ILP-QSI).

QML-IL employs a new approach to QML. It demonstrates that using the incremental decomposition strategy, the learning efficiency can be improved dramatically. This makes QML-IL suitable for learning models of complex dynamic systems, such as biological systems.

A major limitation of QML-IL is that the set of variables and the number of constraints must be provided to each stage by the decomposition solution, and this is difficult in some cases. Although some preliminary work has been carried out on automatic decomposition, further research and experiments are needed to investigate the adaptability and scalability of the decomposition algorithm. In addition, the incremental learner uses the same search strategy as ILP-QSI, so it has the same limitations. This suggests an interesting line of research.

3.10.1 Model representation

Morven and JMorven: Morven (Coghill, Reference Coghill1996) is a state-of-the-art QR engine, which possesses all the benefits of QSIM and introduces many useful features. The introduction of differential planes (Wiegand, Reference Wiegand1991) and vector envisionment (Morgan, Reference Morgan1988; Coghill, Reference Coghill1992) makes it convenient to reason about more than two derivatives. By introducing Fuzzy Set theory (Zadeh, Reference Zadeh1965), each qualitative variable in Morven is associated with a fuzzy quantity space, which makes it possible for Morven to perform fuzzy vector envisionment (Coghill, Reference Coghill1996) and deal with fuzzy qualitative data. JMorven (Bruce & Coghill, Reference Bruce and Coghill2005; Bruce, Reference Bruce2007) is a Java implementation and extension of the original framework of Morven.

Models represented in JMorven: The QSIM and JMorven Model (up to two diffrential planes) for a two-compartment model (described in Figure 2) is shown in Figure 3. In the figure, c1 and c2 denote the concentrations in compartments 1 and 2, respectively; f12 stands for the flow from compartments 1–2; u and f2o are the input and output flows, respectively. In a JMorven model, Plane 0 contains the constraints that construct the actual dynamic system. The constraints in Plane 1 are obtained by differentiating the corresponding constraints in Plane 0. The quantity spaces of the qualitative variables are composed of trapezoidal fuzzy numbers, which are represented by four tuples (see Shen & Leitch (Reference Shen and Leitch1993) or Coghill (Reference Coghill1996) for details of this presentation).

Figure 2

A cascaded two-compartment model

Figure 3

The QSIM and JMorven descriptions of the cascaded two compartment model

3.10.2 Pre-processing phase

The learning process can be divided into two phases: the pre-processing and the backtracking search phases. The pre-processing phase is introduced first.

The pre-processing phase of QML-BKFC includes four steps: (i) Constraint generation: This is similar to GENMODEL, except that given the maximum number of hidden variables, it also generates the qualitative constraints for the hidden variables. (ii) Constraint filtering: All the generated constraints in step (i) are checked for consistency with the input qualitative data. (iii) Precalculation: In this phase, the dependency and conflict relations of the remaining constraints are calculated and recorded. (iv) Constraint set partition: The filtered constraints are divided into several subsets, each of these subsets contains all the defining constraintsFootnote 6 for the same system variableFootnote 7. Then the learning task becomes selecting only one constraint from each of the above subsets to construct a candidate model, then check the validity of this model. The correctness of this selection is guaranteed by a theorem (Pang & Coghill, Reference Pang and Coghill2007b), which states that under the assumptions being made in ILP-QSI (Coghill et al., Reference Coghill, Garrett and King2004, Reference Coghill, Srinivasan and King2008), a well-posed model as defined in ILP-QSI contains one and only one defining constraint for each system variable. By this theorem, the search space can be reduced from exponential to polynomial size.

3.10.3 Learning strategy

After the constraint set partition in the pre-processing phase, as described in section 3.10.2, all the remaining constraints construct the refined search space FCS (Filtered Constraint Set), and these constraints are divided into several subsets, each of which contains all the defining constraints for the same variable. Let DS be a set that takes each of these subsets as an element, denoted by DS = {S ₁,S ₂,…,S_N}, where N is the number of variables (including hidden variables), and for any two elements in DS, |S_i|≤|S_j|, if i ≤ j (i,j ∈ 1…N).

The learning task also becomes that of selecting only one constraint from a different S_i in DS to construct a candidate model, then checking the validity of this model. The backtracking algorithm will incrementally add constraints from different subsets in DS. The subset that has the minimum number of legal constraints will be explored first. The legal constraints are defined as the qualitative constraints in FCS, which do not conflict with all the constraints in the current partial model. Once a new qualitative constraint C_i from a subset, say S_i in DS, is added into the current partial model PM, (PM + C_i) will be checked for validity by the well-posed model constraints; if (PM + C_i) is a full model and passes the well-posed model check, it will be further simulated by JMorven (Bruce, Reference Bruce2007) to test the coverage. If (PM + C_i) fails to pass the check, C_i will be excluded in order to avoid searching for hopeless nodes (the models which include PM + C_i). The algorithm then backtracks to the previous node (PM), and selects the next legal constraint in the same subset S_i to continue the search. If no further legal constraints can be selected in S_i, the algorithm will further backtrack to the previous node, (PM − C_j), where C_j is the most recent constraint that was added into PM.

An auxiliary forward checking method is also performed, that is, when a new constraint is added, all the constraints in FCS that conflict with this new constraint will be ignored in the later search process, as they are illegal constraints for the current partial model. The algorithm will not terminate until the search space has been explored exhaustively.

3.10.4 Advanced experiments

On the basis of the proposed algorithm, some advanced experiments, which extend the experiments in ILP-QSI, were performed to investigate the learning reliability and kernel subsets (as defined in ILP-QSI) under different experimental conditions. These conditions are categorized as (i) the specification of the state variables: Footnote 8fully, partially, or not specifying the state variables; (ii) different number of hidden variables: from zero to the number that leads to unsuccessful learning; (iii) different types of hidden variables, making state or non-state variables become hidden.

Two kinds of two-compartment model were selected for the test: The first model is a closed system with no input and output, but with bidirectional flows; the second model is the cascaded compartmental system as shown in Figure 2. QML-BKFC performed the experiments on these two models under the above-mentioned conditions. The learning reliability and kernel set were obtained for each condition in each model. Some meaningful conclusions are obtained according to the experimental result: (i) The state variables are very important for learning. If some of them become hidden variables, the learning task becomes difficult. (ii) Partially or not specifying the state variables will result in a large search space and may lead to unsuccessful experiments. (iii) The non-state variables have relatively less influence on learning, compared with the state variables. (iv) Too many hidden variables can lead to failure to learn.

3.10.5 Evaluation

QML-BKFC performs a systematic search in the model search space, and is a complete algorithm. It is suitable for small- and medium-size problems. Using this complete algorithm, some reliable experimental results can be obtained and this helps us better understand how different factors can influence the QML. Furthermore, as QML-BKFC uses the Morven framework, it has the potential to deal with fuzzy data.

The major limitation of QML-BKFC is its scalability, compared with QML-CSA to be described in section 4.2. The systematic search will become time consuming or even impractical for large-size problems.

4.1.1 Learning strategy

QME employs a modified GA to search in the model space. The output of QME is a set of QSIM-like models. In order to reduce the complexity, the QDE models under study are simplified as follows: (i) An additional ‘RANGE’ constraint is added to the QSIM constraint repertoire. This constraint is to indicate the legal range of the qualitative variables. (ii) Like GOLEM, the corresponding values are ignored. (iii) In the U-Tube model, the variable that denotes the difference between two levels, together with its related monotonically increasing constraint are also ignored. The simplified U-Tube model is shown in Figure 4. In addition, the following conditions are also assumed: (i) no hidden variables; and (ii) all the landmark values are known beforehand.

Figure 4

The simplified U-Tube model

Encoding strategy: QME uses binary trees to encode the individuals (chromosomes): The leaf nodes denote the qualitative constraints and the branching points represent the hierarchical structure. This encoding method is quite similar to the approach taken in genetic programming (Koza, Reference Koza1992).

Crossover and mutation: The crossover operation is performed on two parents: randomly choosing the crossover points for both parents and exchanging the subtree below the crossover point. The mutation is implemented by randomly choosing a mutation node and replacing the subtree below this node by a randomly generated one.

Fitness calculation: The fitness value for each chromosome is calculated by the following three steps: (i) Raw fitness calculation: Like GOLEM, QME uses both positive examples and hand-generated negative examples. The raw fitness function is the sum of the following two ratios: the ratio of the covered positive examples to all the positive examples, and the ratio of the excluded negative examples to all the given negative examples. If the above two ratios are both equal to 1, an additional bonus score, which is a decreasing function of the model size, will be added to the raw fitness. (ii) Fitness sharing: In order to prevent the premature convergence of the population, raw fitness is shared by similar chromosomes. The phenotype similarity function is used for calculating the shared fitness, that is, the similarity function for any two models is determined by the number of examples (positive and negative) that both models can cover. (iii) Linear scaling: After the fitness sharing, all the fitness values undergo a linear transformation, such that after scaling, the average of these fitness values remains unchanged and the most fit individual is assigned the value of a certain multiple (namely the scaling factor) of the average fitness.

On the basis of the above-mentioned operators, the GA will iteratively explore and exploit the search space in a beam search manner, until the target models are found or the maximum number of generations is exceeded.

4.1.2 Experiments

It is reported that to learn the U-tube, 17 positive and 78 negative examples were used, and QME found five equivalent solutions. Other physical systems, such as the spring system and resistor-capacitor (RC) circuit system, were also tested. In these experiments, QME always found the model that covered all the positive data and excluded all the given negative data.

4.1.3 Evaluation

QME is an evolutionary approach that utilizes a GA to fulfil the learning task. The experimental results show that GA, as a general-purpose searching and learning method, has potential utility in the field of QML. On the basis of the building block hypothesis (Goldberg, Reference Goldberg1989), some meaningful sets of qualitative constraints (building blocks) will emerge if appropriate genetic operations are performed, which makes the search efficient. As the chromosomes vary in size, QME does not need to specify the size of the target model in advance, and in the search process, the algorithm will find the solution model with the appropriate size. The inherent parallelism of GA makes the search more efficient.

Although it benefits from the advantages of GA, QME also suffers from the limitations of GA: (i) It is not a complete algorithm, this means that QME cannot always guarantee to find the target model. (ii) Premature convergence of the population may happen if the related parameters, especially the crossover and mutation probability, are not chosen appropriately. In addition to the influences of GA, there are also other limitations: Like GOLEM, negative examples have to be used; the model has to be simplified by, for example, ignoring the corresponding values.

4.2.1 Evolutionary search

A CSA (de Castro & Von Zuben, Reference de Castro and Von Zuben2000, Reference de Castro and Von Zuben2002; de Castro & Timmis, Reference de Castro and Timmis2002) is used for searching in the model space. CSA is a population-based evolutionary algorithm inspired by the clonal selection theory (Burnet, Reference Burnet1959) of the immune system.

When negative examples are not available and the positive data are complete, the landscape of the model space may have numerous local optima, each of them satisfying the well-posed model constraints, but not all of them generating exactly the same data as the given ones (this means some of the models are underconstrained and some are overconstrained). In addition, when only incomplete data are provided, there may be more than one candidate solution that can cover the incomplete input data. CSA is particularly suitable for dealing with the above situations, because in the above cases, the search space is highly multi-modal and CSA has proven to be efficient for searching in the multi-modal problem space.

The basic idea is that each antibody is encoded to represent a potential model, and in the process of the evolutionary search the antibodies are expected to find the target models.

Encoding strategy: Rather than the binary string encoding used in many conventional evolutionary algorithms, QML-CSA uses an integer encoding strategy. In QML-CSA, each antibody is composed of several slots, and each slot corresponds to a subset, which contains all the defining constraints for the same variable. The value of each slot, which is an integer, indicates a qualitative constraint selected from the corresponding subset. The correctness of this encoding strategy is also guaranteed by the same theorem as was used in QML-BKFC (see section 3.10). In addition, to avoid the antibody getting trapped in a local optimum for an intolerable length of time, a survival time is associated with each antibody. If an antibody exceeds its survival time, it will be replaced by a randomly generated antibody. Figure 5 shows an example of antibody encoding for the two-compartment model in Figure 2. In this antibody, S ₀, S ₁, and S_i contain all the defining constraints for f12, f2o, and the hidden variable Hid0, respectively.

Figure 5

Antibody encoding

Each antibody in the antibody repertoire (population) will be cloned to create several copies (clonal expansion). All the cloned antibodies, forming a temporary population, will be submitted to the ‘hyper-mutation’ scheme. The hyper-mutation is correspondingly modified due to the modification of the antibody encoding strategy. The modified hyper-mutation operator will replace the value of each slot in the antibody with a randomly generated integer, and the range of this randomly generated integer is from 1 to N (N is the number of constraints in the corresponding subset). This means that each constraint in the model will be replaced by a randomly selected constraint within the same defining subset. In the search space, each mutated antibody can be considered as a neighbourhood of the original antibody.

The affinity evaluation will consider the following aspects (as in ILP-QSI): model conflict, dimensional consistency, model language, model connection, model completeness, model singularity, and causal ordering and coverage. These criteria establish a scoring system for evaluating the affinities of the antibodies, such that the model satisfying more criteria get a higher score.

The re-selection process is based on the above scoring system. The antibodies, which have higher scores and have not exceeded their survival time, will be selected to enter the next generation’s population (antibody repertoire).

The modules described above, clonal expansion, hyper-mutation, affinity evaluation, and re-selection, are the basic components of CSA. CSA will iteratively execute these components until it finds the target model or exceeds the maximum number of generations.

4.2.2 Experimental design and results

Two kinds of experiments were designed to test the reliability and scalability of QML-CSA.

• • Reliability tests: For simple models with small sets of training data, the same two kinds of two-compartment models as those used in QML-BKFC were selected. The learning reliability of the algorithm was tested by the power set experiments, as in ILP-QSI. Thus, the kernel subsets and learning curve, defined by Coghill et al. (Reference Coghill, Srinivasan and King2008; see section 3.5.4), were obtained.

• • Scalability tests: For more complex models, a three-compartment and a four-compartment model were selected. They are the extensions of the model in Figure 2. The scalability of QML-CSA was tested by providing complete data but with different hidden variables.

Both of the above experiments were performed by QML-CSA and QML-BKFCFootnote 9. QML-BKFC served as a verification algorithm to test the correctness of QML-CSA for small-size problems, and were also used to compare the efficiency of QML-CSA. The experimental results showed the following two conclusions: (i) For the same power set experiments on small-size models, both QML-BKFC and QML-CSA obtained the same results as in ILP-QSI, that is, the same learning reliability and kernel subsets. (ii) For more complex models, QML-CSA was much more efficient than QML-BKFC, especially when the number of hidden variables was increased. For example, for learning the three-compartment model with three hidden variables, it took QML-BKFC more than 7 days to find the target model, and for QML-CSA, it was approximately 11 hours, which is a 15-fold improvementFootnote 10.

4.2.3 Evaluation

QML-CSA is another evolutionary approach to learn QDE models. Unlike QME, it does not need negative examples, as the dimensional information and well-posed model constraints serve the function of negative examples to some degree. QML-CSA can perform the same kernel set experiments as ILP-QSI and obtain consistent results that demonstrates its reliability. More importantly, compared with QML-BKFC, a major advantage of QML-CSA is its scalability, which enables it to learn large-scale models more efficiently. In addition, like QML-BKFC, it also has the potential ability to deal with fuzzy data because of the adoption of the Morven framework.

The limitations of QML-CSA are as follows: (i) Like ILP-QSI, as it is based on the well-posed model constraints, it cannot learn a model which does not satisfy all the properties of well-posed models; for example, it does not deal with hybrid dynamic systems. (ii) Like QME, QML-CSA is not a complete algorithm; it does not guarantee to find the target model if one exists.

6.3.1 Plant water-balance system

A plant water-balance system was proposed for study by Ramachandran et al. (Reference Ramachandran, Mooney and Kuipers1994). This model describes how the plant balances the amount of water in its body when the amount of water in its external environment changes. Four QDE models were used to model this system. In order to learn this model, two qualitative behaviours, one of which contains 16 qualitative states, were provided to MISQ-RT.

MISQ-RT initially generated eight QDE models. After the unification module, six QDE models were obtained. A further analysis showed that these six models were equivalent to the four QDE models that were generated by the domain experts.

6.3.2 Predator–prey model

Learning the predator–prey model, which is also called the Lotka–Volterra model (Lotka, Reference Lotka1925), was studied by QML-IL (Srinivasan & King, Reference Srinivasan and King2008). The ODEs of the simplified model are as follows:

(12)

(13)

In the above, N and P are the population of prey and predators, respectively; g(N) is the growth rate of the prey, which is a linear function with the constant coefficient; c(P) is the rate of predation upon the prey, and d(P) is the death rate of the predators. If c(P) and d(P) are assumed to be a monotonically increasing function of P, the above ODE model can be converted to a QDE model. Then the qualitative data provided to QML-IL are generated by simulating this QDE model. QML-IL used a one-stage decomposition (which meant not decomposing the model at all) to learn the model. After running QML-IL, five models were returned, one of which was the target model.

6.4.1 Learning the glycolysis pathway

The QDE model of the glycolysis pathway was learnt by ILP-QSI (King et al., Reference King, Garrett and Coghill2005) and QML-IL (Srinivasan & King, Reference Srinivasan and King2008). The glycolysis pathway is one of the most important pathways in biology, which involves 15 metabolites, they are pyruvate (Pyr), glucose (Glc), phosphoenolpyruvate (PEP), fructose 6-phosphate (F6P), glucose 6-phosphate (G6P), dihydroxyacetone phosphate (DHAP), 3-phosphoglycerate (3PG), 1,3-bisphos phoglycerate (1,3BP), fructose 1,6-biphosphate (F16BP), 2-phosphoglycerate (2PG), glyceraldehyde 3-phosphate (G3P), ADP, ATP, NAD, and NADH. As H⁺, H₂O, and orthophosphate are assumed to be ubiquitous in the cell, they are not included in the pathway. Figure 6 shows all the reactions involved in the pathway. The first reaction is taken as an example to illustrate the meaning of these reactions; it means that Glc can be converted to G6P, this reaction is catalysed by the enzyme hexokinase. In this reaction, the ATP will be consumed and ADP will be produced.

Figure 6

The reactions in the glycolysis pathway (from King et al. (Reference King, Garrett and Coghill2005))

Modeling metabolites and enzymes: To model this pathway qualitatively, the enzymes and the metabolites are modelled as components. An enzyme is assumed to have at most two substrates and at most two products. If there are two substrates in a reaction, they are considered to form a substrate complex before the reaction, and the amount of the substrate complex is proportional to the amount of the two substrates multiplied together. It is similar with the products. Each enzyme is considered as a component: the ‘flow’ through this enzyme equals the amount of the substrate complex that enters into it minus the amount of the product complex that leaves from it. Each metabolite is also considered as a component correspondingly. Therefore, the ODE and QDE models for the metabolites and the enzymes are as follows:

Metabolites:

(14)

DERIV(Metabolites, Mdt),

SUM(Flow ₁, …, Flow_n, Mdt).

Enzymes:

(15)

PROD(Metabolite ₁, …, Metabolite _s, S),

PROD(Metabolite ₁, …, Metabolite_p, P),

M⁺(S, Ds),

M⁺(P, Dp),

SUB(Ds, Dp, Flow),

MINUS(Flow, Flow_minus),

where S represents the input substrates and P the products. PROD and SUM are the extensions of MULT and ADD in QSIM, respectively. Based on the above qualitative constraints, more general constraints can be generated to represent the metabolite and enzyme components, as shown below:

where in the ENZYME components, S1 and S2, denote the substrates, P1 and P2 are the products, enzymeFlow is the return value of the component, indicating the flow through this enzyme, which has been described before. While in the METABOLITE component, metaboliteConc and metaboliteFlow are the concentration and flow of this metabolite, respectively, and the third argument denotes all the enzyme flows that can influence this metabolite. The value of the sign can be ‘+’ or ‘−’, which means increase and decrease of this metabolite, respectively. For instance, [enzyme ₁Flow,+] means the flow through enzyme ₁ will increase this metabolite.

Modeling the glycolysis pathway: These two kinds of general constraints (components) are finally used to model the glycolysis pathway, as shown in Figure 7. Two constraints are taken as examples to illustrate the meaning of this model. The first constraint, ENZYME((Glcc, ATPc),(G6Pc,ADPc),Enz1f), means that the flow through Enz1 (hexokinase), Enz1f, controls the transformation of the concentrations of Glc(Glcc) and ATP(ATPc) into G6Pc and ADPc. The 17th constraint, METABOLITE(PEPc,PEPf,(Enz9f, +),(Enz10f, −)), means that the concentration and flow of PEP, PEPc and PEPf, are controlled by the flow through Enz9 (Enolase), Enz9f, and the flow through Enz10 (Pyruvate kinase), Enz10f. Enz9 will produce PEP (according to the sign ‘+’) and Enz10 will remove PEP (according to the sign ‘−’).

Figure 7

The qualitative representation of the glycolysis pathway (from (King et al. (Reference King, Garrett and Coghill2005))

Learning process: First, all the possible ways of combining all the metabolites are generated. Then the domain knowledge of chemical and biochemical reactions is used to reduce the search space. Finally, 18 reactions, which are both chemically and biochemically plausible are selected to construct the initial search space. The qualitative data for learning the pathway are obtained by simulating the target model.

After several months of compute time on a cluster computer, 35 possible models (out of more than 27 000 models) of the glycolysis pathway, among which the target model is included, were returned by ILP-QSI. The long compute time is acknowledged, as it is an approach that would not be done in practice; the problem would generally be tackled in stages. It is a proof-of-concept, and shows that ILP-QSI can narrow down the set of the candidate models to a small scope even when some domain-specific knowledge is omitted.

For QML-IL, the glycolysis pathway was reconstructed in three stages, namely priming, splitting,and phosphorylation. The priming stage includes the first three reactions shown in Figure 6; the splitting stage includes the fourth and fifth reactions; and the rest of the reactions belong to the phosphorylation stage. By this decomposition, QML-IL learnt the model stage by stage. It was reported that the computation times for different stages were 6 seconds in stage 1135 seconds in stage 2 and 5296 seconds in stage 3. After the three-stage learning, only two models were returned, one of which is the target model.

6.4.2 Learning the detoxification pathway of methylglyoxal

QML-Morven was used to learn the qualitative model of the detoxification pathway of methyglyoxal (Pang & Coghill, Reference Pang and Coghill2009). Methylglyoxal is a naturally occurring toxic electrophile, which is harmful to cells. Its detoxification pathway was initially studied by Ferguson et al. (Reference Ferguson, Totemeyer, MacLean and Booth1998), and current understanding of this pathway is shown in Figure 8. In this pathway, the first reaction follows the mass action law, and the second and third are enzymatic reactions, which are taken to obey Michaelis–Menten kinetics (Michaelis & Menten, Reference Michaelis and Menten1913). Glyoxalase I and II are the enzymes that catalyse the second and third reactions, respectively.

Figure 8

The detoxification pathway of methylglyoxal (based on Pang & Coghill (Reference Pang and Coghill2009))

Research on this pathway is still ongoing, and current understanding is based on incomplete knowledge and limited experimental results. Consequently, the quantitative analysis (de Almeida et al., Reference de Almeida, Ozyamak, Miller, de Moura, Booth and Grebogi2008) based on current understanding cannot precisely describe the dynamics of the system. This makes it helpful to reconstruct this pathway at a qualitative level using QML-Morven.

The converting algorithm: On the basis of the standard assumption that all non-enzymatic reactions obey the law of mass action and all enzymatic reactions follow Michaelis–Menten kinetics, a converting algorithm was proposed to convert a possible pathway comprised of several biochemical reactions to QDE models.

Learning strategy: Given known species and types of reactions, all possible reactions are generated. A possible pathway is a subset of all these generated reactions, and all possible pathways constitute the search space. A search algorithm will list the candidate pathways and convert these pathways to QDE models by the converting algorithm. Then the converted QDE models are simulated by JMorven, and the obtained simulation results are compared with the given qualitative data.

Learning results: Given the complete set of 33 qualitative states obtained from simulating the target model, QML-Morven significantly narrowed down the number of possible pathways from 18 772 to 13. The corresponding QDE models for these 13 pathways could generate exactly the same set of states, which was the same as the given complete set.

To deal with the potential complexity of this pathway, a CSA was used to search the target pathway. This CSA is different from the one used in QML-CSA, because the search space for this problem is composed of pathways, rather than qualitative models in QML-CSA. To perform search in this kind of search space, a scoring system was set up to evaluate a pathway, and this scoring system served as a heuristic function to guide the search. Three criteria were used in this scoring system: completeness (the pathway includes all species involved), consistency (no conflicting relations in the pathway), and coverage (the corresponding qualitative model of the pathway covers all the given data).

The CSA was tested not only on the original methylglyoxal detoxification pathway, but also on two more complicated pathways created artificially based on the original one. Experimental results showed that CSA achieved better performance than the exhaustive and backtracking search, which indicated a promising way of improving the scalability of QML-Morven, when learning such complicated systems (Pang & Coghill, Reference Pang and Coghill2009).