IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO.

4, AUGUST 2005 361

Nonlinear System Identification Using

Coevolution of Models and Tests
Josh C. Bongard and Hod Lipson

Abstract—We present a coevolutionary algorithm for inferring

the topology and parameters of a wide range of hidden nonlinear
systems with a minimum of experimentation on the target system.
The algorithm synthesizes an explicit model directly from the
observed data produced by intelligently generated tests. The
algorithm is composed of two coevolving populations. One popu-
lation evolves candidate models that estimate the structure of the
hidden system. The second population evolves informative tests
that either extract new information from the hidden system or
elicit desirable behavior from it. The fitness of candidate models
is their ability to explain behavior of the target system observed in
response to all tests carried out so far; the fitness of candidate tests
is their ability to make the models disagree in their predictions. We
demonstrate the generality of this estimation-exploration algorithm
by applying it to four different problems—grammar induction,
gene network inference, evolutionary robotics, and robot damage
recovery—and discuss how it overcomes several of the pathologies
commonly found in other coevolutionary algorithms. We show
that the algorithm is able to successfully infer and/or manipulate
highly nonlinear hidden systems using very few tests, and that the
benefit of this approach increases as the hidden systems possess
more degrees of freedom, or become more biased or unobservable.
The algorithm provides a systematic method for posing synthesis
or analysis tasks to a coevolutionary system.
Index Terms—Coevolution, evolutionary robotics, gene network Fig. 1. Two approaches to system identification. (a) A batch approach: a set of
inference, grammar induction, nonlinear topological system iden- input data is generated before identification begins; the hidden system generates
tification. output for each input; and the resulting corpus of input/output pairs is used to
construct a model of the hidden system. (b) An iterative approach, in which a
continual model of the hidden system is refined based on each new input/output
I. INTRODUCTION data pair.

S YSTEM identification is a ubiquitous tool in both science

and engineering. Given some partially hidden system, the
experimenter applies a series of intelligently formulated experi-
The batch approach [shown in Fig. 1(a)] involves first gener-
ating a set of input vectors and obtaining the corresponding set
ments to the system in order to learn more about it [44] or make of output vectors, and then generating a model that correctly pre-
it produce a desired output [27]. The experiment is often de- dicts the observed outputs for the given inputs. Many data-inten-
scribed as a set of input to the system, and the resulting be- sive machine-learning methods operate in this way; it is mostly
havior of the system given this input is described by a set of suitable when data is freely available or can easily be collected,
output values. The result of the system identification process is but cannot be controlled. Examples include the evolution of
a model that describes the salient features of the system itself. models to explain stock market data, or models of gene coac-
There are two ways to approach the model inference process: tivation from microarray data.
A batch (offline) approach and an iterative (online) approach. In many application areas, however, data is critically limited,
expensive, or may be biased in unknown ways. In such cases, it
is worthwhile to intelligently formulate new experiments to de-
Manuscript received June 16, 2004; revised November 9, 2004. This
work was supported in part by the U.S. Department of Energy under Grant liberately collect useful data. In the iterative approach [shown in
DE-FG02-01ER45902 and in part by the NASA Program for Research in Fig. 1(b)], the inference process generates an input vector, ob-
Intelligent Systems under Grant NNA04CL10A. This research was conducted tains the resulting output vector, and uses the input/output vector
using the resources of the Cornell Theory Center, which receives funding
from Cornell University, New York State, federal agencies, foundations, and pair to improve the original model. This process is iterated until
corporate partners. a sufficiently accurate model is obtained. Importantly, in the it-
The authors are with the Computational Synthesis Laboratory, Sibley School erative approach the current model can be used to guide the se-
of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY
14853 USA (e-mail: JB382@[Link]; HL274@[Link]). lection of new input vectors that will extract the most new infor-
Digital Object Identifier 10.1109/TEVC.2005.850293 mation about the hidden system given what is already known.
1089-778X/$20.00 © 2005 IEEE
362 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

Fig. 2. Applying the algorithm to different problem domains. (a) The general organization of the estimation-exploration algorithm. (b) Applying the algorithm
to the problem of grammar induction. (c) Applying the algorithm to the problem of gene network inference. (d) Applying the algorithm to evolutionary robotics.
(e) Applying the algorithm for automated recovery after unanticipated robot malfunction.

The field of identification for control [27] involves the gener- and models and tests can be evolved not just to infer a system
ation of a possibly incomplete model solely for the purpose of but also to improve its performance. Here, we describe a number
control, i.e., the correctness of the model itself is important only of implementations of the algorithm and its application to four
insofar as it enables the derivation of a useful controller. How- different problems: inferring finite-state automata; inferring ge-
ever, in this field, as in system identification, data is passively netic regulatory networks; improving behavior transferal from a
collected from the target system; the model is not used to de- simulated to a target robot; and allowing robots to automatically
termine which experiment to perform next. Also, it has been diagnose and recover from unanticipated malfunctions.
noted [27] that no methods exist yet for the automated creation The next section introduces the estimation-exploration algo-
of a controller, given a set of approximate models of the target rithm. The following four sections describe each application in
system. This paper presents an advance toward such automation. turn. The last section provides discusses general properties of
By intelligently choosing input vectors, the number of exper- the algorithm and avenues for further study.
iments that need to be performed on the target system in order
to create a sufficiently accurate model of it can be greatly re- II. ESTIMATION-EXPLORATION ALGORITHM
duced. Similarly, any biases present in batch-generated data can
be actively compensated for by deliberately asking for new tests. The estimation-exploration algorithm is comprised of two
This is important in domains in which it is expensive, risky or populations: The estimation population, which evolves im-
time-consuming to perform experiments, or in which experi- provements to models of the hidden system, given pairs of
ments alter the structure of the hidden system, thus complicating input/output data obtained from the system; and the exploration
the inference process. population, which evolves intelligent tests to perform on the
The method described here is a form of active learning. This hidden target system using the best models so far. A cyclical
machine-learning approach actively seeks out tests that will im- implementation of the algorithm comprises two phases: The
prove the generalization ability of a classifier or learner [16]. estimation phase and the exploration phase. The algorithmic
However, in typical active learning methods, the “model” of the flow of the algorithm is given in Fig. 3.
hidden system is simply a mapping that best translates the input The estimation phase begins with an initial population of can-
data into the observed output data: It does not mirror the internal didate models, which can be random, blank, or seeded with
structure of the hidden system. The algorithm presented here some prior knowledge about the target system. The exploration
synthesizes an explicit model of the hidden system using intel- phase evolves tests that will serve one or possibly two functions:
ligently selected tests. This explicit model can then be used to First, tests should extract as much information about the internal
learn about possible causalities in the hidden system, to localize structure of the system as possible, given what is already known
some change in the system (such as a broken part), or to use that about the system by the candidate models. Second, tests should
localized information to generate desired behavior (such as re- elicit some useful behavior from the target system. Once such
covery of function by circumventing the broken part). a test is evolved, it is applied to the target hidden system and
In this paper, we introduce a coevolutionary algorithm, which some output is obtained. Using this new input/output data pair,
we call the estimation-exploration algorithm, that automates plus any additional input/output pairs obtained during previous
both model inference and the generation of useful experiments. tests, the estimation phase evolves a better model that can better
We use an iterative (online) approach, and maintain two coe- reproduce the observed output data, given the input data. The
volving populations: One that evolves candidate models, and new candidate models are then passed to the exploration phase.
one that evolves tests. The fitness of candidate models is their The cycle continues until either a sufficiently accurate model of
ability to explain the behavior of the target system observed the hidden system is obtained, or a test eventually elicits some
in response to all tests carried out so far; the fitness of candi- desirable behavior from the target system. A schematic repre-
date tests is their ability to make the models disagree in their sentation of the algorithm is shown in Fig. 2(a).
predictions.
The estimation-exploration algorithm can be implemented in A. Coevolution
a variety of ways: The two populations can coevolve in parallel There has been much interest in coevolutionary algorithms
(in steady state) or in a two-phase cycle; models and tests can be within the evolutionary computation community, as evidenced
evaluated against current populations or against their histories, by the recent literature (e.g., [12], [33], [49], [52], [54], and
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 363

Fig. 3. Estimation-exploration algorithm outline.

[58]), starting with the seminal work of Hillis [29] on sorting Pathology I: The Red Queen Effect: One related set of
networks. Contrary to conventional evolutionary systems, co- pathologies derives from the purely subjective measure of
evolutionary systems consist of one or more populations, where fitness in pure coevolutionary systems, and is known as the
individuals may influence the relative ranking of others indi- Red Queen effect: Two populations continuously adapt to each
viduals [12]. For example, whether individual is inferior or other and their subjective fitness improves, but they fail to make
superior to individual may depend on a third individual any consistent progress along the objective metric. Conversely,
rather than on some external fitness metric that provides an ab- they do make progress along an objective fitness but their
solute ranking. There are a number of different forms of co- subjective fitness does not reflect this and falsely indicates lack
evolution: Antagonistic coevolution (e.g., predator–prey), coop- of progress. Both of these effects may initially occur in our
erative coevolution (e.g., symbiosis) or nonsymmetric systems system: A series of difficult tests may decrease the subjective
(e.g., host–parasite or teacher–learner [22]). fitness of models while in fact they are improving in their
In the discussion of coevolutionary systems, it is important to objective fitness. Alternatively, a series of biased tests may give
distinguish between the notion of objective fitness versus sub- rise to overspecialized models which seem to be doing well
jective fitness. Objective fitness is the well-defined absolute fit- at explaining the available data (high subjective fitness) but in
ness metric used in classical evolutionary algorithms. Subjective fact are departing from the true model (lowering their objective
fitness is the fitness as defined by coevolving individuals, which fitness).
may be only weakly correlated with the objective fitness and Both of these effects are transient in our method because of
may sometimes even be misleading. A coevolving individual the anchoring provided by the fixed target system. As the true
only knows its subjective fitness. In the examples presented in objective correctness of models improves, they are ultimately
this paper, we show absolute fitness only for benchmarking pur- able to describe test data and, therefore, their subjective fit-
poses: the algorithm itself has no access to the absolute fitness, ness also increases. Conversely, arbitrary overspecialization is
as in a realistic application, we do not know how close the model removed because it is a source of disagreement among models
really is to the hidden target system; rather it only has indirect and is, therefore, challenged by new tests.
evidence by comparing input–output sets. Pathology II: Cycling and Transitive Dominance: Because
Implementations of coevolution are notoriously difficult, and the subjective fitness criteria is changing over time, individuals
have been plagued with a number of pathologies arising from may “forget” previous abilities, and then rediscover them later,
complex coevolutionary dynamics [34], [62]. Much of the focus yielding a cycling performance. Similarly, an individual may be
of current research has been to address these drawbacks [15], superior to a second individual who is at the same time supe-
[33], [54], [58]. However, hybrid coevolutionary algorithms, rior to it according to a different subjective metric. Cycling and
such as the one proposed in this paper, are less affected by transitive dominance problems are eliminated by the fixed target
these pathologies due to the anchoring effect of the stable target system, because all models need to explain all data so far in ad-
system. dition to any new data.
364 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

Pathology III: Disengagement: Disengagement occurs when Defining representations for inputs and outputs is usually a
one population is entirely superior to another population. The simple matter, as they are typically vectors or matrices which
subjective fitnesses of both populations then become constant indicate the values of variables to be fed into the system or ob-
and all selection pressure is lost, resulting in drift. In system tained as output. Similarity metrics for outputs are often some
identification disengagement may occur when a test may be pre- form of normalized error function, but occasionally more so-
sented which is too difficult for the model population to ex- phisticated metrics may be required as shown in some of our ap-
plain, and all models get an equally low subjective fitness. Al- plications described later. The choices of representation, varia-
though disengagement was not observed in the four applica- tion and comparison of systems usually involve domain-specific
tions presented here, it has occurred in more recent experiments. considerations.
We have recently [6], [63] proposed several mechanisms that 2) Initialization: Initial populations of models and tests need
successfully prevent disengagement from occurring in the esti- to be created. In the absence of any prior information, these
mation-exploration algorithm. The first mechanism is the “test models and tests may be set at random or left blank; if some prior
bank” [6] in which difficult tests are withdrawn from the test knowledge exists, it may be used to bias the initial populations.
suite and are only reintroduced when models become accurate In all four of the applications reported in this paper, the initial
enough to explain them. The second mechanism [63] involves populations are seeded with random models and tests.
searching explicitly for lower difficulty tests by looking for tests 3) Exploration Phase: Useful tests (inputs) are evolved with
which create a less disagreement among models, and reducing their fitness proportional to their ability to create disagreement
the disagreement until the population reengage. Although both among the successful candidate models. Since successful
these methods have been shown to be empirically useful, Ficici models are already compatible with all prior input/output sets,
[64] presents a useful theoretical framework for determining creating disagreement among their predictions focuses the
under which conditions monotonically increasing performance tests on targeting any remaining uncertainties in the model.
can be ensured in any coevolutionary system which attains a This is the approach taken in active learning methods [16].
monotonically increasing knowledge of the search space. Note, however, that disagreement among models can only be
In our algorithm, the evolution of a better model in the esti- measured insofar as candidate models are different; in absence
mation phase allows the exploration phase to evolve a better test. of sufficient model diversity, it may be necessary to seek a
In other words, if little is known about the target system, tests diversity of outputs compared with tests performed on the
must be suggested at random. However, if something is known target system in previous cycles. Either way, we elicit some
about the system, tests can be evolved that cause the system to previously unobserved behavior from the model, and by exten-
exhibit some behavior it has not shown before, generating more sion, from the target system itself.
information about the system. Conversely, the best previously If the overall objective of the entire process is not just to infer
evolved model may fail to replicate a new input/output data pair a system but also to make it behave in a particular way, then the
obtained from the system, thus causing new selection pressure fitness of inputs needs to also capture a measure of its ability
to produce a new model that explains all previous input/output to elicit the desired output. In this case, a form of multiobjec-
data pairs, as well as the new one. tive search may be necessary (e.g., by alternating, weighting,
The power of our algorithm is threefold. First, by evolving or Pareto-selecting tests), though often these objectives coin-
intelligent tests, it is possible to reduce the amount of testing cide. Pareto optimization has already been investigated in the
required on the target system. Second, the algorithm is problem context of coevolutionary algorithms [22], and may be useful
domain independent: the outline of the algorithm given above in evolving tests that extract desirable behavior and information
does not presuppose any particularities about the hidden system about internal structure from the target system.
or the type of experiment to be performed on it. Third, the We have found that in order to facilitate evolution of new
algorithm produces an explicit model of the hidden system. tests, it is often useful to evolve tests from scratch, rather
Fig. 2(b)–(e) sketches the application of the algorithm to the than seed them with tests from the previous cycles. Once
four problems described in the next four sections. successful tests have been found, the best test is carried out on
the target system and the output measured. The input/output set
is recorded with all previous input/output sets.
B. Algorithm Outline
4) Estimation phase: Useful candidate models are evolved
Six steps must be followed to apply the algorithm to a given with their fitness proportional to their compatibility with all
problem, given in Fig. 3. input/output sets collected from the target system so far. As-
1) Characterization of the Target System: This involves suming the target system is consistent, all input/output sets are
defining the target system itself, specifying what aspects of the equally important and should be equally weighted.
target system are known and which aspects must be inferred, When formulating a compatibility error for a specific appli-
and establishing representations for the space of models, inputs, cation of the algorithm, the error often takes the form
and outputs. Variation operators need to be defined to search
the space of models and space of inputs (tests). A similarity
metric comparing two models needs to be established to as-
sess convergence. A similarity metric comparing output also where indicates the number of experiments performed so far
needs to be defined in order to quantify disagreement in model on the target system, is the output obtained from the target
predictions. system during the th test, and is the output obtained from
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 365

TABLE I
OVERVIEW OF THE FOUR APPLICATIONS

the model for the same test. An encoded model that correctly tion phase fails to find any model that is compatible with all ob-
describes the target system will produce the same output as the served data, or when the exploration phase fails to find a test that
target system for all tests, and its error will be low. We call this causes different models to disagree. In the former case, this may
the subjective quality of a model, since it is estimated based indicate that the search space or variation operators do not span
only on known, possibly biased test data; the absolute quality the target system, or that the target system is behaving inconsis-
of a model is not known to the algorithm, but the objective of tently. In the latter case, this may indicate that there is some un-
the exploration phase is to create a suite of tests that make the observable aspect of the target system that the test space cannot
subjective quality approximate the absolute quality. Note that elucidate. In either of these cases, the representation, operators,
diversity maintenance is important since assessment of tests is or similarity metrics need to be reconsidered. In practice, it may
based on creating disagreement among models, and this is only be difficult or impossible to measure the actual accuracy of a
possible insofar as models are different. model: some external kind of validation (see below) may then
In some cases, models may fail to sufficiently explain the be required. However, even if validation is not possible in a prac-
most recent test; in other words, the error of the best model tical situation, the algorithm can still be run for a fixed budget of
at the end of the current pass through the estimation phase will physical trials, depending on the expense of performing a single
be higher than of the best model from the previous pass. This physical trial.
can indicate the beginning of disengagement, which is one of For the grammar induction application, the algorithm is ter-
the three pathologies that coevolutionary algorithms may expe- minated if either criteria 1) or 2) is met: either a perfect model is
rience [15], [33], [54], [58], [62]. Disengagement occurs when discovered, or 10 model evaluations have been performed.
one population poses too much of a challenge to the other pop- For the gene network inference application, criterion 3) is used:
ulation; the dominated population then loses its fitness gradient, the algorithm terminates when 10 model evaluations have
because all individuals perform equally badly against the indi- been performed. For the final two robot applications, the algo-
viduals from the dominating population. In recent papers [6], rithm terminates when a fixed number of target evaluations have
[63], we have proposed several mechanisms for combating dis- been performed [criterion 3)].
engagement. The most useful mechanism has proven to be the 6) Validation: In inference applications, a cross validation
test bank, in which difficult tests are removed from the test suite step helps assess the significance of the resulting successful
and only returned to the test suite when models become accu- model. Cross validation is achieved by performing a previously
rate enough to explain them. However, in the four applications unseen test on the target system and comparing its outputs to
presented here, disengagement was not detected, so these mech- the prediction proposed by the model. If the validation step is
anisms are not used. unsuccessful, the new input/output set is added to all previous
5) Termination: There are four criteria for terminating the input/output sets, and the algorithm resumes at the estimation
algorithm, either 1) a sufficiently accurate model of the target phase.
system has been obtained; 2) an evolved test has caused the These six steps are applied, in turn, to each of the four prob-
target system to exhibit the desired behavior; 3) a maximum lems described in this paper. An overview of the applications is
number of target or model evaluations have been performed; or given in Table I. In the next section, we describe the application
4) the algorithm failed. The algorithm fails when the estima- of the algorithm to the problem of grammar induction.
366 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

III. APPLICATION 1: GRAMMAR INDUCTION

Grammar induction is a special case of the larger problem

domain of inductive learning [5], which aims to construct sys-
tems based on sets of positive and negative classifications (for
an overview, see [50]). In grammar induction, the system is a
finite-state machine, or FSM, that takes strings of symbols as Fig. 4. Example of deterministic finite automaton. The DFA is composed of
three states (n = 3), assumes a binary alphabet, and the third state is a final
input, and produces a binary output, indicating whether that state. The first of the two sentences receives a positive classification because it
string, or sentence, is a part of the FSM’s encoded language. ends at state 2; the second sentence receives a negative classification because it
Grammar induction is an attempt to model the internal struc- ends at state 1.
ture of the hidden FSM, based only on pairs of sentences and
classifications. The most successful grammar induction algo- which is given as . For simplicity, we here choose a binary al-
rithms produced so far are heuristic in nature (see [14] for an phabet. We then can represent as a matrix with values
overview), and there have been several attempts to evolve FSMs in , where each column in represents a state, and the
from sample data [25], [45], [46], but all approaches so far as- values in a given column indicate which new state to transition
sume that a representative sample of sentences have already to, given the current symbol of a sentence. Finally, we can rep-
been classified by the target FSM. In other words they follow the resent the set of final states using a binary string of length :
batch-oriented approach to system identification [see Fig. 1(a)]: if , then state is a final state; otherwise, it is not. For
a pregenerated set of input/output data is fed into some induc- simplicity, we assume that the first column of represents the
tion algorithm, and a model FSM is produced. start state. For the example three-state DFA shown in Fig. 4, the
Moreover, the pregenerated set of sentences in these ap- sentence 010 receives a positive classification, but the second
proaches are usually assumed to be balanced: that is, there is sentence 011 receives a negative classification.
an equal number of positively and negatively classified sen- A simple grammar induction problem is then to model the
tences. However, if an FSM is unbalanced (the FSM classifies hidden DFA by finding a and an such that for all sentences,
a majority of sentences either positively or negatively), then the model produces the same classifications as those produced
collecting a balanced set of classified data for induction may by some target DFA. Therefore, from among ,
be prohibitively expensive: a large number of sentences would , , and are assumed to be known, but and are as-
have to be classified before a sufficient number of sentences sumed to be unknown. Note that a model with and/or
attaining the minority classification could be collected. The ap- may still reproduce all of the same classifications as the
plication of our algorithm to the problem of grammar induction target DFA, if it accurately captures the classification structure
does not presuppose a set of input/output data, but rather uses of the target DFA using a different and , or if only a few
an evolved set of models of the hidden FSM to generate a new input/output pairs have been tested.
sentence which, when fed into the hidden FSM, will produce In the estimation phase, each genome represents a candidate
a new data point that will help improve the model. It will be model of the hidden DFA. Since the number of states is
shown in this section how the estimation-exploration algorithm known, each genome encodes a integer matrix with
dramatically reduces the number of sentence classifications values in and a binary vector of length . If
required by the target FSM, compared with a batch approach to sentences have already been classified by the hidden DFA, then
grammar induction. the quality of a candidate model—its subjective error—is given
A deterministic finite automaton, or DFA, is a type of FSM by
that can be represented using the five-tuple ,
where is the number of states, is the alphabet of the (1)
encoded language, is a transition function, is the start state,
and is a set of final states. Then, given some sentence made where is the binary classification generated for the th
up of symbols taken from the alphabet , and beginning at sentence by the hidden target DFA, and is the classification
the start state , the first symbol is extracted from the sentence, generated by the model DFA for the same sentence. Then, a
and based on that symbol the sentence transitions to a new genome that obtains perfectly reproduces all of
state as indicated by . The next symbol is then extracted from the classifications produced by the hidden DFA so far, and
the sentence, and based on the sentence transitions to a new genomes with higher values of reproduce less of the
state. This process is continued until all symbols in the sentence correct classifications.
have been exhausted. If the last state visited is a member of , The space of inputs (tests) that can be applied to any given
then the sentence receives a positive classification (the sentence DFA is the space of binary sentences. For simplicity, we only
belongs to the language); otherwise, a negative classification is evolve sentences with length , the number of states. Sentences
assigned (the sentence does not belong to the language). Given of variable length could also be used. Thus, each genome in the
this formulation, we can now apply the estimation-exploration exploration phase is a binary string of length .
algorithm to this problem. The quality of a sentence is determined to be the ability of
1) Characterization of the Target System: We assume that the the given sentence to cause the greatest variation in classifica-
target system is a DFA, and that we know the number of states, tions produced by the algorithm’s current best set of candidate
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 367

models. In order to evolve more than one accurate model during When all 50 generations have been evaluated, the sentence
the estimation phase, variation must be maintained in the model with the lowest value of is output and supplied to the target
population until the estimation phase terminates. The simplest system.
way to ensure this is to evolve isolated subpopulations, so that 4) Estimation Phase: The estimation phase, like the ex-
each subpopulation will evolve an accurate and unique model. ploration phase, evolves a population of 200 initially random
So for this application, the population of candidate models is genomes for 50 generations. However, the population of the
partitioned into subpopulations during each pass through the estimation phase is partitioned into two genetically isolated
estimation phase. At the end of each pass, there are then can- subpopulations with 100 genomes each.
didate models that explain the hidden DFA: the best model from Also, like the exploration phase, when all genomes in each
each subpopulation obtained during that pass through the esti- subpopulation have been evaluated, selection and mutation oc-
mation phase. The quality metric of a candidate sentence (ex- curs. In each subpopulation, 150 pairs of genomes are selected
periment) is then given by at random, and for each pair the genome with the lower value of
is copied and replaces the genome with the higher value
of . The copied genome then undergoes mutation: first
(2)
either the encoded or is selected with equal probability;
then a value within the selected matrix or vector is randomly se-
where is the classification of the candidate sentence by lected and replaced with a new random value. Crossover is also
candidate model . Sentences that do not differentiate between not implemented in this phase. The newly created genomes are
the candidate models—all models produce the same classi- evaluated, and evolution continues. When this phase terminates,
fication—obtain (poorest quality); sentences that two competing candidate models are output to the exploration
produce the maximum classification variance obtain phase: the model with the lowest value of from each sub-
(best quality). When a sentence is evolved that induces high population.
classification variance, and that sentence is classified by the Unlike the exploration phase, which begins each pass with
hidden DFA, then the resulting classification will usually a population of random genomes, on the second and subse-
lend support to candidate models during the next pass quent passes through the estimation phase each subpopulation
through the estimation phase, and provide evidence against the of random genomes is seeded with the best genome it evolved
remaining half. This idea is borrowed from the coevolutionary during the previous pass. This allows the algorithm to continue
literature, in which it has been shown that the fitness of a improving its previous best candidate models, and for the results
test should be proportional to its ability to induce a learning of the newly evolved test to lend support or provide evidence
gradient in the competing population of learners [12], [22], against these models.
[33] (in our framework, a learner is a model DFA). The value 5) Termination: After each pass through the estimation phase,
of causing disagreement between models is discussed in more the algorithm is terminated if either of two criteria are met: the
detail in Section VII. most fit model achieves an absolute error of 0, or 10 model
2) Initialization: Since we do not assume an initial model of evaluations have been performed. In this paper, hidden DFAs
a given hidden DFA, the algorithm initially generates a random with ten states were inferred. If a model DFA classi-
binary sentence, which is classified by the hidden DFA. This fies all possible sentences correctly, then it obtains an absolute
single sentence/classification pair is then fed into the estimation error of 0. The set of all possible sentences in this context is con-
phase in order to generate an initial model. sidered to be all binary strings with lengths equal to the number
3) Exploration Phase: The exploration phase begins each pass of states in the target DFA. So for the inference of DFAs with
with a population of 200 randomly generated genomes, and ten states, there are 2 possible sentences. The number
the population evolves for 50 generations. At the end of each of model evaluations is tallied as followed: during the th pass
generation, when all of the genomes have been evaluated, 150 through the estimation phase, each model requires evaluations
pairs of genomes are selected randomly and sequentially (some in order to determine its fitness; during the exploration phase,
genomes may be selected more than once), and for each pair the each test requires model evaluations in order to determine its
genome with a lower value of is copied and replaces the fitness, where is the number of model subpopulations.
genome with a higher value of . 6) Validation: The calculation of the absolute error of a model
The copied genome is then mutated: one randomly selected DFA is considered to be validation, because it uses sentences
element from the candidate sentence is chosen with a uni- unseen during model evolution. Thus, from the viewpoint of
form distribution,1 and is replaced with a new random value. validation, the algorithm can be said to have been successful
Crossover is currently not used, but may be implemented in if the absolute error of the model DFA output after termination
future improvements to the algorithm. A total of 150 replace- is below some sufficiently small error threshold.
ments are performed after each generation. Note that a given
genome may undergo more than one mutation if it is selected, A. Grammar Induction Results and Discussion
copied, mutated, and then selected again. Once a set of selec- The second column in Table I provides a summary of the
tions, replacements, and mutations have occurred, all of the application of the estimation-exploration algorithm to this par-
new genomes in the population are evaluated. ticular instance of grammar induction. The proposed algorithm
1As are all the other random values used in the experiments described in this was run against a number of target DFAs with ten states, along
paper. with two control algorithms. The first control algorithm was a
368 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

Fig. 6. Mean performance of the three algorithms against a single target DFA.
Forty independent runs were performed using each algorithm against the same
target DFA composed of ten states. Absolute errors were computed for the best
model DFAs produced at the end of each generation; the errors are averaged
over the 40 runs.

DFAs in that algorithm. Fig. 5(b) reports the absolute errors for
the same model DFAs.
As can be seen, only the estimation-exploration algorithm
achieves a perfect model, after about 190 target evaluations; the
other two control algorithms terminate without ever finding a
perfect model. Also, Fig. 5(a) shows that the standard GA be-
gins with a very high subjective error, while the latter two al-
Fig. 5. Model accuracy results from a single run. (a) Subjective errors of the gorithms begin with zero subjective errors, which then climb
most fit model DFAs produced by the two control algorithms (a standard GA and later reapproach zero. This is because it is much easier for
and random testing) and the proposed algorithm (intelligent testing). (b) The the model DFAs to correctly classify the few sentences labeled
absolute errors for the sample model DFAs. The most fit model DFA at the end
of each generation, for all three algorithms, was extracted and its subjective and by the target DFA than later on, when many sentences have
absolute error calculated. been labeled. Fig. 5(b) shows that the standard GA exhibits
the fastest drop in absolute error, which is to be expected be-
standard genetic algorithm: all 1024 binary sentences of length cause early in the algorithm the standard GA has much more
10 were first classified by the target DFA, and then a random observed data available for model evolution than the latter two
population of 200 model DFAs were evolved until either a per- algorithms do. Finally, Fig. 5(a) shows that intelligent testing
fect model is discovered (one that achieves an absolute error induces much greater subjective errors in the model DFAs com-
of 0) or 10 model evaluations are performed. Genetic en- pared with random testing. This is indirect evidence that the pro-
coding, selection and mutation are the same as described in posed algorithm is discovering informative tests: tests that ex-
the previous section. This control algorithm, thus, follows the pose unexplained aspects of the target system.
batch approach to system identification depicted in Fig. 1(a). Fig. 6 reports the mean absolute errors of the most fit model
The second control algorithm is identical to the estimation-ex- DFAs produced by all three algorithms for one of the target
ploration algorithm, except that the exploration phase is dis- DFAs. The absolute errors for the model DFAs were averaged
abled: when the estimation phase terminates, a randomly gen- over the 40 independent runs. The error bars in the figure (as
erated binary sentence of length 10 is sent to the target DFA for well as the error bars in all subsequent plots in this paper) in-
classification. dicate standard error with a 95% confidence interval. As in the
Thirty target DFAs with ten states each were randomly gener- previous plot, the average accuracies of the model DFAs pro-
ated for inference. Each of the three algorithms then performed duced by the standard GA early on are much better than the
40 independent runs against each of the 30 target DFAs, re- latter two algorithms, due to the scarcity of observed data avail-
quiring a total of independent runs. able to the iterative methods. However, the proposed algorithm
Fig. 5(a) reports the subjective errors [see (1)] of the best attains a significantly better mean performance than the other
model DFA during the course of each algorithm’s execution two algorithms after less than 100 target evaluations have been
against the same target DFA. For the standard genetic algorithm performed.
(GA), the subjective error of the most fit model DFA after each Fig. 7 reports the comparative performance of the three al-
generation was recorded; for the latter two algorithms, the sub- gorithms against the 30 target DFAs. Fig. 7(a) reports the mean
jective error of the most fit model DFA after each generation of number of model evaluations required by each algorithm to pro-
the estimation phase terminates was recorded. Note that because duce a perfect model: a model DFA with an absolute error of 0.
the standard GA uses all possible sentences for model evalua- Note that in most cases, the proposed algorithm is able to con-
tion, subjective error is equal to absolute error for the model sistently discover a perfect model with fewer model evaluations
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 369

Fig. 7. Comparison of algorithm performance against 30 random target DFAs. (a) The mean number of model evaluations required until a perfect model is found.
(b) The mean number of target evaluations required until a perfect model is found. (c) The mean fraction of positive classifications acquired by the sentences
proposed to the target DFA.

than the standard GA, even though the standard GA has access DFA 5 shown in Fig. 7(c)]. It can be seen that intelligent testing
to much more information about the target system (compare the tends to consistently elicit more balanced classifications from
heights of the dark gray bars against the light gray bars). Con- very imbalanced target DFAs: in other words intelligent testing
versely, the standard GA always outperforms iterative, random creates a sufficiently accurate model that more of the underrep-
testing (compare the dark gray bars against the medium gray resented class of sentences can be discovered and proposed to
bars), indicating the performance gain of the estimation-explo- the target DFA (indicated by the relatively taller light gray bars
ration algorithm comes not from performing iterative inference, on the left side of Fig. 7 and the relatively shorter light gray bars
but from sending informative tests to the target DFA. on the right side.) This indicates that evolving tests that cause
Fig. 7(b) compares the number of target trials required by the disagreement among models leads to the uncovering of less ob-
second control and proposed algorithms to find a perfect model. servable components of the target system: in this applications,
The standard GA is not shown in this plot because each stan- this involves traversal to final states that produce minority clas-
dard GA requires exactly 1024 target trials before the algorithm sifications.
commences. For most of the 30 target DFAs, intelligent testing However, it can be seen that intelligent testing also outper-
requires significantly fewer target trials in order to discover a forms random testing on many balanced DFAs (note the signif-
perfect model. Even with roughly 1/5 the number of target trials icant performance advantages of the proposed algorithm on the
(about 200 trials compared with 1024 trials), the proposed algo- horizontally centered target DFAs). This indicates that a target
rithm can discover a perfect model more quickly (i.e. with fewer DFA may have less observable components even if it is bal-
model evaluations) than the standard GA can. anced: for example, some states may be visited by random sen-
Fig. 7(c) compares the fraction of positive classifications ob- tences much more often than other states if the state transition
tained from the target DFA for the tests proposed by the second graph has regions of dense connectivity and other regions have
control and proposed algorithms. Because the second control al- sparse connectivity.
gorithm proposes random tests, the fraction of positive classifi- In this section, we have described the application of the esti-
cations obtained by this algorithm gives an approximation of the mation-exploration algorithm to the problem of grammar induc-
balance of the target DFA. The balance of a DFA is considered [Link] have shown that it is worthwhile to evolve experiments
to be the probability that it will produce a positive classification (sentences) for this problem, and that the quality of a test lies in
given any arbitrary sentence; DFAs with probabilities of 0.5 are its ability to distinguish between two candidate models. Further-
considered balanced, and those with probabilities far from 0.5 more, we have shown that intelligent testing can discover a per-
are considered imbalanced. As can be seen, the target DFAs are fect model using fewer model and target tests than either a sim-
aligned in Fig. 7 in order of their balance. [The target DFAs in ilar iterative algorithm that employs random testing, as well as a
Fig. 7(a) and (b) are aligned with the ordering in Fig. 7(c); for standard genetic algorithm that assumes a large amount of pre-
example, target DFA 5 in Fig. 7(a) and (b) is the same target classified sample data. Finally, we have shown that intelligent
370 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

testing is most useful when the hidden DFA has low observ-
ability. In the next section, we describe the application of the
algorithm to the problem of gene network inference.

IV. APPLICATION 2: GENE NETWORK INFERENCE

Systems biology [39] is concerned with the synthesis of large Fig. 8. Observability of gene regulation. In the gene network model used
here, it is assumed that gene product concentration is measured after some time
amounts of biological detail in order to infer the structure of period has elapsed (t ). (a) If the final concentration of some gene product is
complex structures with many interacting parts. In this way, sys- intermediate, there is only one possible rate for the gene product concentration
tems biology can be viewed as another example of system iden- change between t and t (we assume in this paper that the change during this
time interval is linear). (b) If the final which attempts to evolve an accurate
tification. The field of gene network inference is a rapidly bur- model using a set of observed gene product concentration is either zero or
geoning subfield in systems biology [39], and is concerned with completely saturated, there are several hypotheses for how this gene is regulated
inferring genetic regulatory networks based on the results of a (three possible hypotheses are indicated by trajectories i, ii, and iii). Therefore,
the regulation of the gene in (b) is less observable than the gene in (a).
set of tests performed on the network in question.
Many different models of the underlying genetic network
have been used, usually classified based on the amount of Both the target gene network and models of it are represented
biological detail inherent in the model (see [17] and [18] for an using an matrix with entries in [ 1, 1]. The new
overview). In addition to the type of model, several methods concentration of gene after some time period is then given by
have been used to infer genetic networks, including clustering
algorithms (see [18] for an overview), correlation metrics (3)
[3], linear algebra [13], simulated annealing [48], and genetic
algorithms [31], [38].
A number of input and output data pairs are required in order Since the variables indicate concentration, the and
to obtain enough information about the target network to infer functions bound the value between 0 (for no concentration)
its structure correctly. As pointed out in [18], it is desirable to and 1 (for concentration saturation). Genomes in the estimation
reduce the number of input/output pairs required as much as phase are then represented as matrices. The subjective
possible, so that a minimum of experiments have to be con- error of a given network model is simply the mean squared
ducted. Also, the type of experiment required should be as cheap error between the concentrations output by the hidden network
as possible in terms of experimental difficulty and accuracy of and the candidate model so far. For this application, absolute
acquired output data. Iba and Mimura [31] showed that by using error is simply the distance of a model gene network from the
a multipopulation evolutionary algorithm to not only infer the actual gene network
hidden network, but also to propose additional experiments that
would most help to refine the current best evolved network hy- (4)
pothesis. However, their model requires the experimenter to per-
form costly knockout or lesion experiments in order to supply where represents gene ’s regulation of gene in the target
the algorithm with an actual subset of the regulatory network. system , and represents gene ’s regulation of gene in the
Our approach does not require such invasive experiments, but candidate model .
rather assumes that the input data to a biological experiment is a The test for a hidden gene network is determined as a vector
set of chemical concentrations. These chemicals can be viewed of floating-point values in [0, 1] that represent initial gene
as either initial concentrations of the gene products themselves, product concentrations; thus, genomes in the exploration phase
or some other media such as signalling proteins or glucose that are vectors of floating-point values. The output from the target
externally affect the cell. This approach obviates the need for in- system or model is a vector of floating-point values indicating
vasive testing, leading to a simpler experiment. The input is as- eventual gene product concentrations after some time period has
sumed to trigger gene regulation, which produces a set of output elapsed. Unlike in the grammar induction application in which
gene product concentrations over time, as obtained using mi- test quality is given by its ability to cause the candidate models
croarray tools. We can now apply the estimation-exploration al- to disagree, the test quality in this application is to increase the
gorithm to infer the structure of a hidden gene network, given observability of the system; in other words, to minimize the
sets of initial chemical concentrations, and resulting product number of extremal concentrations, which indicate less about
concentrations, as outlined in Fig. 2(c). the internal structure of the network than intermediate concen-
1) Characterization of the Target System: The target system trations. Fig. 8 shows why final gene product concentrations that
is assumed to be a gene network in which the number of genes are either zero or saturated reveal less information about the un-
is known, but how one gene contributes to the regulation of derlying network than gene product concentrations that fall be-
another is assumed to be unknown. tween these extremes.
The steps for applying the estimation-exploration algorithm 2) Initialization: The algorithm begins by generating some
to this particular problem are only summarized below: please initial random input vector . This is then applied against the
refer to [9] for more details regarding this application. hidden target network . The resulting input/output
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 371

vector pair is fed into the estimation phase, and the algorithm only gets through 55 iterations before the 10 model eval-
commences. uations have been exhausted). The standard GA was provided
3) Exploration Phase: The exploration phase begins with with 100 sets of sample data so that it requires almost twice as
a population of 1000 random genomes, and the population is many target evaluations as the other two algorithms.
evolved against the single candidate model output by the es- All three algorithms were run against a total of 80 randomly
timation phase. Selection and mutation then proceed similarly generated target gene networks. The networks varied in the
to the previous application, with the mutation operators altered number of genes in the network , and the network’s connec-
to handle the floating-point valued genomes. The population is tivity . For each network, each gene was regulated by a total
evolved for 30 generations, at which point the genome with the of other randomly chosen genes (with the possibility of itself
best quality is output to the target network . being one of the genes); the other genes do not directly
4) Estimation Phase: Evolutionary search in the estimation regulate that gene . The first 20 target networks were
phase begins with a population of 200 random genomes, each generated with and ; the following 20 with
encoding some . The subjective error of each genome is and ; the following 20 with and ; and last
computed, and 150 replacements are performed: for each 20 networks with and . For each target network,
genome pair, the genome with lower subjective error replaces each algorithm performed a total of 30 independent runs in an
the genome with higher subjective error. Mutation but not attempt to infer the hidden regulatory network.
crossover is applied to the copied genomes. The population is Fig. 9 reports the average performance of the three algorithms
evolved for 30 generations, at which point the genome with against these 80 target networks. Because a perfect model is
the lowest subjective error is passed to the exploration phase. never found, each algorithm performs 10 model evalua-
During the second and subsequent passes through this phase, tions; the standard GA performs 100 target trials; the iterative
the initial random population is seeded with the best model algorithm with random testing performs 56 target trials; and
evolved during the previous pass. the estimation-exploration algorithm performs 55 target trials.
5) Termination: As in the grammar induction, the algorithm The mean performance of each algorithm is calculated to be the
was run until either a perfect model was produced, or a max- mean of the absolute errors of the 30 best model gene networks
imum number of model evaluations (for this application 10 ) output at the termination of the 30 independent runs.
had been performed. Because the model in this application con- First, it can be noted that for all classes of target network, the
sists of continuous values (indicated by the strength and kind of standard GA performs significantly worse than the other two al-
gene ’s regulation of gene ), there is a vanishingly small prob- gorithms. Individual runs (data not shown) indicate that the stan-
ability that a perfect model—one that achieves an absolute error dard GA tends to converge on very poor local minima after only
of 0—will be discovered, so all runs for this application termi- a few generations, and never make any subsequent progress.
nate when 10 model evaluations have been performed. Surprisingly, the iterative algorithm employing random testing
6) Validation: Validation of a candidate model is interpreted performs significantly better than the standard GA, using half
as its absolute error. as many target trials. This indicates that the iterative approach
for this problem domain provides an evolutionary advantage:
A. Gene Network Inference Results and Discussion each successive pass through the estimation phase presents the
In order to test the algorithm against this task, two control evolving models with a fitness landscape different from the pre-
algorithms were formulated as in the previous application: the vious pass, and that this difference is induced by the new test
first control algorithm acts as a standard GA which attempts to obtained from the target system.
evolve an accurate model using a set of observed data from the Second, the proposed algorithm only outperforms random
target system; the second control algorithm is identical to the testing on the largest and most densely interconnected gene net-
estimation-exploration algorithm, except that the exploration works ( , ). This occurs for two reasons. Because
phase is disabled: each pass through the exploration simply the tests are composed of continuous rather than binary or in-
outputs a random set of gene product concentrations. teger values, there is a very low probability that the same test
Because the tests in this application are composed of con- will be proposed more than once during random testing, and a
tinuous values (instead of binary values in the grammar induc- set of unique tests is more informative than a set of tests in which
tion application), there are an infinite number of possible tests, several tests appear more than once. Also, large, densely inter-
so we can not provide the standard GA with an exhaustive set connected networks tend to lead to more extremal gene product
of sample data already processed by the target system. Rather, concentrations than smaller, sparse gene networks: a gene reg-
we generate a set of 100 random vectors of initial gene ulated by many other genes will exhibit a large change in its
product concentrations, with the values of each vector sampled gene product concentration, compared with a gene which is only
from [0, 1]. By choosing a maximal number of model evalu- regulated by a few genes. Thus, as in the previous application,
ations 10 , and using a population size of 200 for all our algorithm is most valuable for inferring large target systems
three algorithms, and running both the estimation and explo- with low observability.
ration phases for 30 generations, the second control algorithms This claim that extremal gene product concentrations are
performs a total of 56 target evaluations, and the proposed al- more likely in large, dense gene networks is supported by
gorithm performs 55 target evaluations (because the proposed Fig. 10: 1000 random networks were generated using values
algorithm incurs slightly more model evaluations per iteration of selected from [2, 30] with a uniform distribution, and
than the second algorithm due to the active exploration phase, it values of selected from , where is an already
372 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

Fig. 10. Observability of various gene network types. 1000 gene networks
with different values of n and k were generated randomly. For each network,
100 random input vectors were supplied, and the average fraction of resulting
extremal output concentrations was calculated. (a) The relationship between the
number of genes (n) and the fraction of extremal values. (b) The relationship
between connection density (n=k ) and the fraction of extremal values for the
same 1000 networks.

Fig. 11. Testing results from a sample run using random and intelligent
testing. Each column indicates the resulting gene product concentrations
obtained from a test proposed by either (a) the iterative algorithm using random
testing or (b) the proposed algorithm. Black squares indicate extremal gene
product concentrations and white areas indicate nonextremal concentrations.

the corresponding ten output vectors were calculated using

(3). The average fraction of extremal output concentrations
was computed for each network, and is plotted in Fig. 10(a)
as a function of the number of genes in the network , and
in Fig. 10(b) as a function of connectivity . Clearly, the
fraction of noninformative output elements increases both with
the number of genes, and with connectivity. Thus, it becomes
increasingly valuable not only to evolve models, but also to
Fig. 9. Mean performances of the three algorithms against the four classes of evolve informative experiments for the inference of larger and
target gene networks. (a) Mean performance of all three algorithms against the
20 target gene networks with n = 5 genes and k = 2 connectivity, (b) the 20 more complex gene networks.
target networks with n = 5 and k = 5, (c) the 20 target networks with n = 10 Fig. 11 substantiates this claim by showing the resulting gene
and k = 2, and (d) the 20 target networks with n = 10 and k = 10. product concentrations from tests proposed by a sample run of
the iterative algorithm with random testing and the proposed al-
randomly selected value for . For each of the 1000 random gorithm. The sample runs were both drawn from the set per-
networks, ten input vectors were randomly constructed, and formed against one of the target gene networks with and
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 373

Fig. 12. Mean fraction of extremal gene product concentrations produced by

random and intelligent testing. Results are reported for both algorithms against Fig. 13. Evolutionary progress of a traditional evolutionary robotics regime.
the 20 target gene networks with n = 10 and k = 10. Thirty generations of a genetic algorithm that evolves a neural network
controller for a robot so that it performs a forward locomotion task (the robot
and its controller are shown in Fig. 14). Fitness is the distance traveled by
. As can be seen, the proposed algorithm tends to pro- the robot. Each dot indicates a simulated robot evaluation. The best controller
is then transferred to a target robot, which has some hidden morphological
pose tests that produce many fewer extremal gene product con- difference (one of the robot’s lower legs has broken off). Evolution then
centrations [indicated by the fewer black squares in Fig. 11(b) continues on the target robot, using the same genetic algorithm, for another 70
compared with Fig. 11(a)]. generations. Each square indicates an evaluation on the target robot. The solid
line indicates the best fitness achieved at each generation.
Fig. 12 reports the mean fraction of extremal gene product
concentrations produced by the tests proposed by the iterative
algorithm with random testing and the estimation-exploration The following two sections describe the application of our
algorithm, while inferring the 20 target networks with algorithm to two important problems in robotics research:
and . It is clear that for all 20 target networks, intel- automating behavior generation, and automating damage diag-
ligent testing achieves significantly fewer extremal concentra- nosis and recovery.
tions than random testing. This indicates that evolving tests that
produce as many nonextremal gene product concentrations as V. APPLICATION 3: EVOLUTIONARY ROBOTICS
possible allows for the evolution of more accurate tests when An evolutionary robotics experiment requires an evolutionary
the target network is large and has low observability. algorithm to optimize aspects of a simulated or physical robot in
Interestingly, some of the genes in a given hidden network order to generate some desired behavior. Because it is difficult
are less observable than others. For example the experimental and time-consuming to perform the thousands of fitness evalua-
results obtained by the proposed algorithm in Fig. 11(b) fail to tions required by an evolutionary algorithm on a physical robot,
consistently achieve nonextremal concentrations for the fourth most or all of evolutionary robotics experiments are performed
gene (indicated by the long black bars for that gene over the in simulation.
course of the run), while it seems relatively easy to obtain nonex- Evolution in simulation raises a major challenge: The trans-
tremal concentrations for the sixth and seventh gene (evidenced feral of evolved controllers from simulated to physical robots,
by the long white bars for those genes near the end of the run). or “crossing the reality gap” [32]. Because there are always dif-
This indicates that this gene is very strongly positively or nega- ferences between the simulated robot and its physical instan-
tively regulated and, thus, has low observability because it often tiation, controllers evolved in simulation do not always allow
reaches extremal concentrations. We plan to improve the quality for the same behavior to arise in the target robot when trans-
metric for experiments in future work such that genes with low ferred. Fig. 13 shows the evolutionary progress for an evolu-
observability (as indicated by previous experiments) automati- tionary robotics experiment. A neural network controller was
cally come under increasing scrutiny. evolved for a simulated legged robot for 30 generations, and
Thus, the estimation-exploration algorithm has three benefits then transferred into a different robot (the target robot, which in
for gene network inference. First, it does not require invasive, this case was also simulated) that is morphologically different
expensive, slow and disruptive experiments such as knockout or from the simulated robot: the target robot is missing one of its
lesion studies. Rather, the exploration phase carefully evolves a lower legs. The transferal causes a large drop in fitness, and
low-cost experiment (a change in the initial gene product con- subsequent evolution, involving about 3000 evaluations on the
centrations) that yields a large amount of information about the target robot, only recovers about 70% functionality.
target system. Second, the number of experiments performed on There are several approaches to the challenge of controller
the target system is reduced, because each proposed experiment transferal, including adding noise to the simulated robot’s sen-
is carefully chosen. Finally, the automated evolution of informa- sors [32]; adding generic safety margins to the simulated objects
tive experiments becomes increasingly valuable as the hidden comprising the physical system [26]; evolving directly on the
networks become larger and more densely interconnected, since physical system ([23], [47] and [59]); evolving first in simula-
large, dense networks often produce information-poor output tion followed by further adaptation on the physical robot ([47],
data in response to random testing. [51]); or implementing some neural plasticity that allows the
374 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

Fig. 14. Robot and its controller used for the evolutionary robotics experiments. (a) Quadrupedal robot. T indicates touch sensors; A indicates angle sensors;
M indicates motorized joints. (b) Neural network controller. Sensors are arranged on the input layer, and motors are arranged on the output layer. B indicates
bias neurons that output a value of 1 at each time step.

physical robot to adapt during its lifetime to novel environments the evolution is performed in simulation or the real world.
([19], [24], [60]). Other approaches have widened evolution’s control over the
Another approach that can be used in lieu of, or in addition design process by subjugating the controller topology and/or
to the above-mentioned approaches, is to treat the problem as a the robot’s morphology to modification as well (e.g., [2],
system identification task in which there are hidden differences [7], [30], [41], [43], and [55]) with the aid of simulation.
between the target robot and the robot simulator, which must In this application, we evolve the simulator itself: this may
be automatically uncovered and included into the simulation. involve virtual modifications to the simulated robot’s body, its
In this section, we document the application of our algorithm sensor/motor apparatus, its virtual environment, or the physical
to this problem, such that the exploration phase evolves con- parameters of the simulation itself. The ability to evolve a sim-
trollers for a sensor-driven target robot to make it behave, and ulated robot’s morphology or its environment, in addition to its
the resulting behavior of the target robot is used to refine the controller, has become much easier recently due to the advent
robot simulator. The goal is to automatically refine the simulator and availability of physics-based simulators, which allow for
sufficiently that controllers evolving in it cause the target robot faster than real-time, three-dimensional dynamic evaluation of
to produce similar behavior to that seen in simulation. Thus, different physical systems (see [21] for an overview). Next, we
the target robot serves as the target system; the controller is the describe the preparatory steps for applying our algorithm to
input vector that elicits behavior from the system; the resulting this problem.
sensor time series are treated as the output from the system; and 1) Characterization of the Target System: The target system
the robot simulator serves as the model of the target system. here is a quadrupedal robot with an articulated body, a set of
Fig. 2(d) outlines the algorithmic flow for this application. sensors and motors, and a neural network controller that con-
Most evolutionary robotics experiments evolve controller nects sensors to motors. In this work, the target robot is not a
parameters for a robot with a fixed controller topology and physical robot out in the world, but a separate, simulated robot
fixed morphology (examples include [23] and [53]), whether which is identical to the default simulated robot except for some
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 375

unknown morphological differences. The task of the algorithm tional information will be extracted from the robot automati-
is to indirectly infer these differences, and modify the default cally. This obviates the need to explicitly include a term in the
simulator to reflect them accurately. In future work, we plan to fitness function to evolve controllers that distinguish between
apply our algorithm to a physical robot. The layout of the target competing models (as in the grammar induction case) or expose
robot is shown in Fig. 14(a), and the topology of its controller is some unobservable part of the target system (as in the gene net-
shown in Fig. 14(b). The simulation environment, the simulated work inference case). However, such terms could be formulated
robot and its neural network controller are described in more and included in the fitness function in future.
detail in [8]. 2) Initialization: Unlike the previous two applications, here
The experimenter must then choose those morphological or the algorithm begins in the exploration phase, evolving neural
environmental characteristics that he thinks may differ between network controllers using a default robot simulator.
the target robot and the simulation and that will affect the suc- 3) Exploration Phase: Each pass through this phase begins
cessful transferal of behavior. Currently, this selection is done ad with a random population of 100 genomes. Each controller in
hoc, but it is intuitive: for example the color between the robots turn is evaluated on the simulated robot, and the fitness (for-
will not matter, but mass distribution and sensor behavior will ward displacement) is recorded. Selection and mutation is im-
most likely not be modeled accurately at the outset, and both plemented differently from the previous two applications.
affect behavior greatly. Once all of the genomes in the population have been evalu-
For the work reported here, we induced a difference in mass ated, they are sorted in order of forward displacement, and the 50
distribution and sensor time lags between the target robot and genomes with the least displacement are deleted from the popu-
the initial default robot simulator. For these initial experiments, lation. Fifty new genomes are selected to replace them from the
we assume that the target robot only differs from the default remaining 50, using tournament selection, with a tournament
simulation in those characteristics that we have chosen to place size of 3. Each floating-point value of a copied genome has a
under evolutionary control. From similar experiments [8], we 1% chance of undergoing a point mutation (replacement of the
have found that in some cases even an approximate simulator evolved value with a new random value). Of the 50 copied and
that does not refine all the differing physical characteristics be- mutated genomes, 12 pairs are randomly selected and undergo
tween the simulated and target robot allows for adequate transfer one-point crossover. The population is evolved for 30 genera-
of behavior. tions. When this phase terminates, the controller with the best
Thus, for this application the algorithm must automatically fitness is output and downloaded to the target robot, and the re-
infer the differences in body part masses and sensor time lags sulting sensor values are recorded. Both the evolved controller
only using sensor feedback from the target robot. Therefore, a and resulting sensor time series are passed into the estimation
candidate model is a set of 17 floating-point parameters (masses phase.
of the nine body parts and the time lags for the eight sensors), 4) Estimation Phase: The genetic algorithm operating in the
which are collected into genomes and evolved. The parameters estimation phase is similar to the one in the exploration phase.
are used to modify the default robot simulator, producing a can- During the first pass through the estimation phase, a population
didate robot simulator. of 100 random genomes is generated; during subsequent passes
The quality of a candidate simulator is given by the ability through this phase, the initial population contains 99 random
of the simulated robot to mimic the observed behavior of the genomes, and a copy of the best genome evolved during the
target robot. More specifically, given previously evolved con- previous pass. Each genome is evaluated in turn: the default
trollers tested on the target robot, the simulated robot should simulation is modified according to the genome; the simulated
mimic as closely as possible the sets of sensor time series pro- robot is evaluated using the previously evolved controllers for
duced by the target robot. However, unlike the previous two ap- 1000 time steps; sets of sensor time series are obtained; and the
plications, we can not directly compare model and target output genome’s rolling mean is calculated. Selection and mutation for
to gauge model accuracy. Quantitatively comparing sensor data the estimation phase is similar to that for the exploration phase.
from two highly coupled, highly nonlinear machines like the The population of the estimation phase is evolved for 30 gen-
robot used here is very difficult: slight differences between the erations as well, and outputs the most accurate simulator to the
two machines rapidly leads to uncorrelated signals. To address exploration phase, which begins again with a random popula-
this, we have formulated a comparison metric called the rolling tion of controllers.
mean metric, which is described in detail in [8]. 5) Termination: The algorithm iterates through the cycle
In this application, a test is a set of 68 floating-point values shown in Fig. 2(d) 20 times, starting at the exploration phase
used to label the robot’s neural network controller. Therefore, with the default simulator, and terminating after the 20th pass
the genomes in the exploration phase are vectors of 68 floating- through the estimation phase. This produces 20 evolved con-
point values in [ 1, 1]. The quality of a test is how far the robot, trollers and 20 simulator modifications.
using the labeled controller, moves forward during 1000 time 6) Validation: There is no validation phase necessary for this
steps of the simulation. It is implicitly assumed that there are example since in this case we know the absolute error. This
many possible gaits for the target robot and, therefore, different error is not available to the algorithm.
evolved controllers that make it move in various ways. There-
fore, as long as separate passes through the exploration phase A. Evolutionary Robotics Results and Discussion
begin with different random populations of controllers, different Fifty independent runs of the algorithm were conducted
useful controllers will be output to the target robot, and addi- against the target robot. Fig. 15 shows the convergence toward
376 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

Fig. 15. Convergence toward the physical characteristics of the target robot. Each pass through the estimation phase produces a set of mass changes for each of
the nine body parts of the robot (top row) and a set of time lags for each of the eight sensors of the robot (bottom row). The trajectories indicate the mean and
standard deviation of the best guesses output by the estimation phase from the 50 independent runs. The each bar corresponds to the best guesses from that pass
through the estimation phase. The open circles indicate the actual differences between the target robot and the starting default simulated robot (for example, the
first body part of the target robot is 1.5 kg heavier than the corresponding body part of the default simulated robot).

Fig. 16. Behavior recovery after controller transferal. After the first pass through the exploration phase, the best evolved controller caused the simulated robot
to move as shown in (a) the trajectory of its center of mass is given by the thin line. The same controller was then supplied to the target robot, and the resulting
trajectory of its motion is given by the thick line. The movement of the updated simulated robot after the 20th pass through the exploration phase (using the new
best evolved controller) is given by the thin line in (b). The motion of the target robot using the same controller is given by the thick line in (b). The horizontal axis
indicates forward distance, and the vertical axis indicates height (both are in meters).

the actual mass distribution and sensor time lags of the target this case, an asymmetric robot) are less observable than other
robot. The figure makes clear that for all 50 runs, the algorithm components.
was better able to infer the time lags of the eight sensors than One of the runs was selected at random, and the gait of the
the mass increases of the nine body parts (indicated by the con- simulated robot was compared against the gait of the target
vergence of the means and the negligible standard deviations robot, when both used the same evolved controller. Fig. 16(a) in-
toward the actual time lags on the bottom row compared with dicates the change in behaviors when the first evolved controller
the large standard deviations on the top row). In the inference was transferred, and Fig. 16(b) shows the behavior change when
of the time lags, the algorithm only had difficulty with sensor the 20th evolved controller was transferred, during the last iter-
3, and tended to overestimate its particular lag. Also, the algo- ation through the algorithm’s cycle.
rithm tended to have difficulty converging on the correct mass
increases of most of the body parts, but did a relatively good This figure shows that even with an approximate description
job inferring the mass changes of body parts 4 and 8. of the robot’s mass distribution, the simulator is improved
The fact that the algorithm had more difficulty with some enough to allow smooth transfer of controllers from simulation
parts of the robot and not others is not surprising: the asym- to the target robot. Using the default, approximate simula-
metric mass distribution of the robot causes some aspects of the tion, there is a complete failure of transferal, as indicated
robot’s morphology to be more observable than others. For ex- by Fig. 16(a): the target robot simply moves randomly, and
ample, the time lag of a touch sensor, attached to a very heavy achieves no appreciable forward locomotion.
leg that is only ever dragged along the ground plane (such that After 20 iterations through the algorithm, an improved sim-
the sensor always fires), is much more difficult to infer than an- ulator is available to the exploration phase, which evolves a
other touch sensor attached to a leg that repeatedly touches the controller that allows the simulated robot to move forward, al-
ground and lifts off again. Also, a sensor with a shorter time lag though not as far as the original simulated robot [indicated by
is probably easier to infer than one that has a long time lag. the shorter trajectory in Fig. 16(b) compared with Fig. 16(a)].
Another reason that it was easier for the algorithm to infer This is most likely due to the fact that the simulator is now ac-
sensor rather than body part characteristics is that the sensors curately modeling the heavier mass of the target robot. Also, the
themselves provide feedback about the robot. In other words, new gait causes the robot to hop [indicated by the large vertical
the algorithm automatically, and after only a few target trials, curves of the robot’s center of mass in Fig. 16(b)] instead of
deduces the correct time lags of the target robot’s sensors, but walk [indicated by the steady trajectory of Fig. 16(a)]. In con-
is less successful at indirectly inferring the masses of the body trast to the first pass, the target robot exhibits very similar be-
parts using the sensor data. As discussed in the previous two havior to the simulated robot when it uses the same controller:
sections, often some of the components in a target system (in both travel a similar distance (about 6.5 m), and both move in
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 377

Some work has been done on employing evolutionary algo-

rithms to restore functionality after some unanticipated damage
has occurred, but all of this work relies on massive numbers of
hardware trials: robot recovery has been demonstrated in [4] and
[47], and for electronic circuits in [10] and [37]. However, re-
peated generate-and-test algorithms for robotics is not desirable
for several reasons: repeated trials may exacerbate damage and
drain limited energy; long periods of time are required for re-
peated hardware trials; damage may require rapid compensation
(e.g., power drain due to coverage of solar panels); and repeated
trials continuously change the state of the robot, making damage
Fig. 17. Average transferal success after each target trial. The light gray bars diagnosis difficult.
indicate the average distance traveled by the simulated robot using the best
evolved controller output by that pass through the exploration phase, over all
Grefenstette and Ramsey [28] propose an algorithm for con-
50 runs. The dark gray bars indicate the average distance traveled by the target tinual behavioral learning in a robot context which does not re-
robot using the same controller, during each target trial. Error bars indicate two quire large numbers of target trials. However, their system as-
units of standard deviation.
sumes that changes in the robot’s environment (in their case,
speed change of a prey robot) can be directly inferred using
the same way (both exhibit a hopping gait that produces trajec- sensor data and provided to the behavior learning component. In
tories with similar frequencies and amplitudes). their case, it is assumed that the prey robot may change speed,
Finally, for each pass through the exploration phase, the dis- and their described algorithm assumes that speed change can be
tance traveled by the simulated robot using the best evolved con- detected by the predator robot; the observed change in then fed
troller was averaged over all the 50 runs. Similarly, the distance directly into the predator’s behavior learning component. How-
traveled by the target robot using the same controller was aver- ever, the more difficult aspect of damage diagnosis and recovery
aged over the 50 runs. The results are shown in Fig. 17, which is usually diagnosis: in the estimation-exploration algorithm it
indicates that improvement in behavior transferal success, even is assumed that robot or environmental changes are not known
using an approximate simulator, is a general phenomenon. On in advance, and must be inferred indirectly using the sensors
average, over the 50 independent runs, there is a drop by 50% that also drive behavior. Grefenstette and Ramsey state that in-
in the distance traveled by the target robot, compared with the direct sensing of environmental change could be included into
default simulated robot. After about five iterations through the their system, but they do not provide any description as to how
algorithm’s cycle there is only a statistically insignificant de- this could be accomplished.
As mentioned in the previous section, several types of plastic
crease in distance traveled between the two robots. Although
neural network controllers have been proposed that allow for
not shown in Fig. 17, this similar distance is matched in all of
rapid, lifetime adaptation to external perturbation (e.g., [19],
the cases viewed by a qualitative similarity in gait patterns, as
[24] and [60]). Furthermore, Keymeulen et al. [36] have formu-
shown for a single run in Fig. 16(b).
lated an algorithm that continuously updates an internal model
In this section, we have shown that the estimation-exploration
of sensor input/world response data obtained from a wheeled
algorithm can be used to automatically evolve an accurate robot
robot, and uses this model to evolve and download controllers
simulator. Even if the simulator does not describe the target
to the robot during task execution. They have demonstrated that
robot perfectly, it does allow for automated transferal of evolved
their algorithm greatly reduces the number of required hardware
controllers from simulation to the target robot. Moreover, suc-
trials, compared with a similar model-free algorithm. However,
cessful transferal becomes possible after only three or four target none of these approaches generate a hypothesis describing what
trials. In the next section, we extend this approach to automate particular change the robot or its environment has experienced.
the process of diagnosis and recovery for robots that undergo Second, the evolved controllers have not been shown capable of
some unanticipated malfunction. fundamental reorganization when faced with unanticipated, cat-
astrophic failure (such as the separation of a limb), as is demon-
VI. APPLICATION 4: DAMAGED ROBOT RECOVERY strated here for our proposed algorithm.
Srinivas [57] was one of the first researchers to study error di-
For a robot to function for long periods of time in a hos- agnosis and recovery, but his approach, along with subsequent
tile, unknown or remote environment, it must be able to deal approaches ([1], [11], [20], [35], [61]), require online operation
autonomously with uncertainty: specifically, unanticipated in- (repeated testing on the physical robot), and cannot handle unan-
ternal damage or external environmental change. The recent dif- ticipated errors. Baydar and Saitou [4] proposed the first offline
ficulties with Jet Propulsion Laboratory’s two Mars rovers pro- error diagnostic and recovery system, which relies on Bayesian
vide a dramatic example: both robots suffered different, unantic- inference for error diagnosis, and genetic programming [40] for
ipated partial failures [42]. Automatic recovery is most acute in error recovery. However, their algorithm also only handles pre-
such instances where human operators cannot manually repair specified error types.
or provide compensation for failure. In this work, we are con- Mahdavi and Bentley [47] recently demonstrated an online
cerned with catastrophic, highly nonlinear robot faults that re- evolutionary algorithm that automatically recovers behavior for
quire recovery controllers qualitatively different from the orig- a physical robot. However, after damage the physical robot re-
inal controller. quired 400 hardware trials and nearly seven hours to recover
378 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

Fig. 18. Parsing of an example genome in the estimation phase. This particular genome has one active and three inactive genes. The resulting single modification
(a weakening of the third motor by 43%) is applied to the default simulator.

72% of its original functionality. Their approach is similar to modification of the simulated robot and its environment; this
that portrayed in Fig. 13: evolution is continued on the target simulation is treated as the candidate model. The quality of a
robot after transferal. model is again the ability of the given simulated robot to mimic
Here, we describe the application of the estimation-explo- the behavior of the crippled target robot.
ration algorithm to the problem of damaged robot recovery, This encoding gives a search space of 200
which, for the results presented here, requires only ten target 10 different genomes. As the values are discretized,
trials (on average) to restore functionality. Like the previous there are at most
application, the algorithm automatically evolves two separate
structures: a robot simulator that explains unanticipated internal possible sets of simulator modifications, although some of
damage or external environmental change suffered by the target these are not unique because ordering is irrelevant, multiple
robot; and a compensatory neural network controller that re- modifications can be additive (two weakenings of the same
stores functionality to the target robot, given the evolved robot motor by 10% is equivalent to one weakening of the same
simulator. Below, we outline the application of the first five motor by 20%) and modifications encoded by inactive genes
steps of the proposed algorithm to the problem of automated are not applied.
damage diagnosis and recovery; validation, as in the previous A test in this application is the same as for the previous ap-
applications, is not used here. plication: genomes encode a vector of 68 synaptic weights for
1) Characterization of the Target System: Like the previous the controller. The quality of a labeled controller is how far it
application, the target robot is an independent simulated robot enables the crippled robot to move forward during 1000 time
that possesses some unknown differences in contrast to an steps of the simulation.
initial, approximate simulator. We use the same quadrupedal 2) Initialization: As in the previous application, we begin with
robot as before (Fig. 14), but the target robot no longer has a a default simulation of the quadrupedal robot and its environ-
larger mass and time-delayed sensors, but does suffer some un- ment, and initiate the exploration phase to evolve a neural net-
expected malfunction or environmental change. The algorithm work controller for it. The best evolved controller is then down-
must then diagnose this malfunction using sensory feedback as loaded onto the damaged target robot, and the resulting sensory
before, and update the simulator to reflect this damage: if the data, along with the evolved controller, are fed into the estima-
target robot has lost a leg, then a good simulator must simulate tion phase in order to improve the simulation.
a robot with the correct leg missing. 3) Exploration Phase: During each pass through the explo-
The model in this application is a genome that can encode ration phase a population of 200 genomes is evolved for 40
zero to four simultaneous malfunctions. Each encoded malfunc- generations. Before each evaluation, the default simulator is
tion is parameterized as to: which body part, joint, sensor, or modified (i.e., the simulated robot is damaged) according to the
motor it applies; what has happened (the joint has broken or evolved set of modifications supplied by the estimation phase.
the sensor has failed); and how severe the damage is (expressed Selection and mutation are implemented the same way as in the
as a percentage). Genomes are comprised of 16 floating-point previous application. At the end of this phase, the best controller
values, and each genome is divided into four genes, each com- is output to the target robot.
prised of four values; each gene describes a possible malfunc- 4) Estimation Phase: Each pass through the estimation phase
tion. These values are parsed into a set of malfunctions to apply begins with a random population of 200 genomes; on the second
to the default simulated robot. Fig. 18 outlines a sample genome and subsequent passes through this pass, the random population
being parsed; refer to [8] for more details regarding the specifics is seeded with the best simulator modification evolved during
of the genetic encoding. Thus, each genome produces some the previous pass. Each genome in turn is parsed; the resulting
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 379

TABLE II
UNANTICIPATED SITUATIONS TESTED

modifications (if any) are applied to the default simulator; the

modified simulated robot is evaluated using the controllers al-
ready tested on the target robot; and both the sensor logs from
the damaged target robot and those obtained from the simula-
tion are compared using the rolling mean metric [8] to assess
the accuracy of the current genome. Once all of the genomes in Fig. 19. (a) The evolutionary progress of a sample run against situation 5:
the population have been evaluated, the population is sorted ac- one of the angle sensors breaks by 50%. The dotted lines indicate the progress
cording to accuracy, and selection, mutation and crossover is ap- of the first five passes through exploration phase. The captions indicate the best
simulator modifications evolved by four passes through the estimation phase.
plied. Each pass through the estimation phase evolves the pop- The triangle shows the original distance traveled by the target robot. The circle
ulation for 40 generations. At the end of the pass, the best set of indicates distance traveled after suffering the failure. The squares indicate
simulator modifications is passed to the exploration phase. the distance traveled by the damaged target robot during each of the four
trials. (b)–(d) Images of the target robot’s motion: (b) before the unanticipated
5) Termination: Starting at the exploration phase, the algo- situation; (c) after encountering the unanticipated situation; and (d) during the
rithm iterates through the cycle shown in Fig. 2(e) ten times, fifth target trial.
leading to ten evolved controllers, ten target tests, and ten
evolved simulations. See Table I for a summary of the applica-
tion of the estimation-exploration algorithm to this problem.
6) Validation: There is no validation phase necessary for this
example since in this case we know the absolute error. This
error is not available to the algorithm.

A. Damaged Robot Recovery Results and Discussion

Thirteen different target robots were tested: each suffered a
different unanticipated situation. Table II describes the situa-
tions. The algorithm was then run 20 times against each of the
13 target robots. The estimation phase attempts to evolve simu-
lator modifications that will describe this situation: for each of
the 13 situations, there are only a few genomes (or none, in the Fig. 20. Automated recovery after a catastrophic failure using the proposed
algorithm. (a) The behavior of the default simulated robot using the best
case of situations 11–13) that will describe it perfectly, out of
controller evolved during the first pass through the exploration phase. (b) The
the 86 400 possible situations. Each run requires about 4 h of target robot has lost part of a leg through joint separation, and the robot attempts
computation on a 1 GHz PC. to move using the first evolved controller. (c) The motion of the damaged target
Fig. 19 documents the performance of a single run of the robot during the second target trial.
control algorithm against a target robot that has suffered a 50%
failure of one of its touch sensors (situation 5 in Table II). The the second test, the algorithm is able to evolve a compensatory
algorithm evolves an approximate description of the damage controller for the damaged robot that restores functionality.
after the first target test, and evolves a perfect description after Fig. 20 indicates that the proposed algorithm is even able to
the second test. This indicates that as the algorithm accumulates provide compensatory controllers when the target robot suffers
more data from the target robot by causing it to behave in dif- catastrophic failure, such as a change in the robot’s topology
ferent ways, it is better able to model the target robot’s situation. (one of the lower legs has broken off, represented by situation
Also, Fig. 19(a) shows that, armed with an accurate model after 3 from Table II). Here, the algorithm correctly identifies the
380 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

accurate robot simulator, and automated diagnosis and recovery

from unanticipated robot malfunction. These four applications
have not only demonstrated that our approach is problem do-
main independent, but have also highlighted its various proper-
ties and benefits. Before we summarize the results from these
applications, however, it is worthwhile to go into more detail
about the differences between the proposed algorithm and other
coevolutionary algorithms.

A. Differences From Pure Coevolutionary Algorithms

Fig. 21. The average performance of the proposed algorithm against all 13 The main difference between the estimation-exploration al-
damaged robots. The algorithm was run ten times against each situation. The
white bars indicate the average distance traveled by the default simulated robot gorithm and other coevolutionary systems proposed so far is that
using the first evolved controller. The light gray bars indicate the average our algorithm is a hybrid system: the fitness of both the model
distance traveled by the target robot during the first target trial. The dark gray and test populations is influenced by the opposing population,
bars indicate the average distance traveled by the target robot during the tenth
target trial. but also by the behavior of an additional, third component—the
target system. In other coevolutionary algorithms, the fitness of
an individual is influenced by one or more individuals from an-
damage after the second pass through the estimation phase, and other (or the same) population. For this reason, the fitness of
the exploration phase uses this accurate simulator to evolve a a coevolving individual in such systems is purely subjective; it
new controller that allows the damaged robot to move using only is calculated based only on other individuals. However, in our
three legs. Note the difference in trajectories shown in Fig. 20(a) system, there is a static component to search, which is the ability
and (c): the original gait caused the robot to move using long of a model to mimic the observed behavior of a target system:
hops [indicated by the long loops in Fig. 20(a)], but the damaged as long as the target system does not alter its internal structure
robot uses a compensatory controller that causes it to move in over time, there is an objective aspect to the search for better
short, hesitant hops] indicated by the short loops in Fig. 20(b)]. models. For this reason, we view previous coevolutionary sys-
Thus, the algorithm automatically diagnoses mild or severe mal- tems as pure coevolutionary systems, and the estimation-explo-
function, and either modifies the gait slightly over subsequent ration algorithm as a hybrid coevolutionary system.
trials or evolves a qualitatively different gait, depending on the
damage incurred. B. Application of the Algorithm
Fig. 21 reports the ability of the proposed algorithm to recover
In the grammar induction problem, it was shown that an in-
function for damaged robots. The algorithm was able to restore
telligent test is often one that distinguishes between alterna-
functionality (regarded as a statistically significant increase in
tive candidate models, a concept related to that propounded in
mean performance between the first and tenth hardware trial) for
the competitive coevolutionary literature: the quality of an indi-
all but damage scenarios 6, 11, and 13. Even in these three cases,
vidual in one population is proportional to its ability to induce a
many of the runs produce a significant increase in performance.
performance gradient in the subset of individuals challenging it
Note that even in the catastrophic cases such as 3, 4, 7, and 10
from the competing population [12], [22], [33]. There are sev-
(indicated by almost complete lack of forward motion during the
eral reasons cited as to why such a test is desirable, but there is
first hardware trial) some functionality was restored, indicating
that our algorithm may be useful for recovering some function a specific advantage in the context of the estimation-exploration
from severely damaged robots. algorithm, described next.
Note also that both algorithms were able to successfully re- If a test does not induce disagreement between models, and
cover when faced with an unanticipated environmental change: the test is passed to the target DFA for classification, then on
the horizontal canting of the ground plane (situation 10). This the subsequent pass through the estimation phase there will be
was probably due to the drastic effect on all sensors induced by an absolute but not a relative change in the fitnesses of the best
this change. It is important to note that the proposed algorithm models: the fitnesses of the best models from each subpopu-
can distinguish between internal damage and external environ- lation will either all go up by one (if they match the classifica-
mental change, and reconstruct the change in simulation, based tion of the target DFA) or down by one (if they disagree with
solely on indirect sensory information. the target DFA). If the best models all improve in fitness, there
Most encouragingly, some functionality was restored in cases is little chance for evolutionary improvement: new models can
11 and 12, in which there was no genome that could perfectly only classify any previous sentences not yet correctly classified
describe the malfunction. Further study is required in order to by the current best models, if there are any. If the best models all
understand how truly unanticipated situations can best be ap- experience a drop in fitness, then there is sufficient room for im-
proximated by a set of simulator modifications. provement: a new model may correctly classify all the sentences
that the best model(s) do, as well as the new sentence that they
do not.
VII. DISCUSSION
However, as the inference process proceeds, and the models
The previous four sections have described the application of in the estimation phase improve, there is an increasing likeli-
the estimation-exploration algorithm to four separate problems: hood that the best models will correctly classify any incoming
grammar induction, gene network inference, the evolution of an sentence. Therefore, there will be a diminishing return over the
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 381

inference process of proposing sentences for classification that and our approach could be combined with other approaches
the best models already agree on. However, if a sentence is pro- that have already been suggested for automating the process of
posed for classification by the target DFA that causes disagree- evolved controller transferal.
ment among the best models, there is a guarantee that of the In the last application, we employed the estimation-explo-
best models will incorrectly classify the new sentence, thereby ration algorithm to automatically diagnose unanticipated mal-
ensuring the possibility of further model improvement. functions suffered by the target robot, and automatically evolve
Tests that cause model disagreement have an added benefit: compensatory controllers for it. As in the previous three appli-
it helps to uncover less observable components of the target cations, recovery was achieved with a small number of target
system. For example, in the grammar induction application, in- trials. In this application a maximum of ten target trials were suf-
telligent testing tended to elicit much more balanced classifi- ficient to restore functionality in the face of: minor, catastrophic,
cations from unbalanced target DFAs than random testing was or compound damages; no damage at all; changes in the robot’s
able to [see Fig. 7(c)]. For many of these unbalanced DFAs, the environment; or malfunctions that could not be perfectly ex-
proposed algorithm was able to extract enough of the minority pressed in simulation. Most importantly, the algorithm was able
classifications using intelligent testing to evolve a much more to automatically distinguish and describe these diverse situa-
accurate model, on average, than the same algorithm but with tions using only sensory feedback and the rolling mean metric,
random testing. In such cases when the hidden system has low which quantitatively compares behaviors between highly cou-
observability, tests must be formulated that produce more of the pled, highly nonlinear dynamic robots. Once a correct diag-
underrepresented output in order to infer the internal structure nosis is discovered the algorithm is often able to evolve a com-
of the system: in the estimation-exploration algorithm this oc- pensatory controller, often evolving qualitatively different gaits
curs as a natural result of evolving tests that cause disagreement in response to fundamental changes experienced by the target
among models. robot.
For the gene network inference problem, we have shown that
it is possible to automatically evolve the regulatory connec- C. Model Accuracy Versus Structural Accuracy
tions between a set of genes given a set of input and output When performing system identification, it is possible to gen-
gene product concentrations. Moreover, our approach does not erate a model that successfully mimics the input–output be-
require invasive, expensive, slow, and disruptive experiments havior of the target system, which we refer to as model accu-
such as knockout or lesion studies. Rather, the exploration phase racy, but fails to reflect the actual internal structure of the target
carefully evolves a low-cost experiment (a set of initial chemical system, which we call structural accuracy. Whether our pro-
concentrations that trigger gene regulation) that yields a large posed algorithm actually captures structural as well as model
amount of information about the physical system. accuracy depends on the specific application.
Like the grammar induction problem, it was found that For grammar induction, there are often several DFAs that pro-
evolving useful tests, instead of proposing random ones, speeds duce identical classifications. So for this application, it is impos-
the discovery of accurate gene network models, thus requiring sible to determine, using any algorithm, whether a model DFA
a minimum of physical experimentation. Moreover, it was has captured the actual internal structure of the target DFA: ab-
found that the usefulness of evolving tests increases as problem solute error of a candidate model is simply approximated using
difficulty increases: it was found that evolving tests is increas- a large set of previously unseen sentences.
ingly useful as the hidden gene networks increase in size and In the gene network application, our algorithm does capture
connectivity. This finding bodes well for the scalability of our the internal structure of the target gene networks: absolute error
approach, suggesting that our algorithm may prove very useful is considered to be the distance between the internal structures
for genetic network inference in particular, and for systems of the model and target networks. As can be seen in Fig. 9(d),
biology research in general. the internal structures of the best model gene networks produced
In Section V, we described the application of the estima- by the estimation-exploration algorithm are closer to those of
tion-exploration algorithm to the problem of automatically the target networks than those models produced by the other
generating controllers for robots. In that application the hidden two algorithms (indicated by the final absolute errors of the best
target system is a robot with unknown morphological differ- models). However, this correlation between model and struc-
ences compared with an initial default simulated robot. The tural accuracy depends on the way in which we have chosen
algorithm evolves controllers for the robot in simulation, which to model gene networks: validation methods for determining
are then tested on the target robot. The resulting sensor data, whether a given gene network model actually captures the regu-
along with the evolved controller, serve as input/output data lation inherent in an actual biological gene network is currently
pairs which are used to evolve improvements to the robot sim- being investigated.
ulator. It was shown that this approach automatically produces In the case of the first robot application, Fig. 15 indicates that
accurate enough simulators for the successful transferal of con- all of the runs tend to converge on the actual internal structure
trollers from simulation to the target robot, using a minimum of of the target system: the mass distribution and sensor time lags
trials on the target robot. This approach is appealing because: of the target robot. In the case of the second robot application,
repeated testing on a physical robot is time-consuming and Table II reports that for many of the damage scenarios, the esti-
potentially damaging; the sensors already being used to drive mation-exploration algorithm actually discovers the perfect de-
behavior also act as diagnostic sensors, returning indirect infor- scription of the damage (see column 3). However, for these two
mation about the target robot’s morphological characteristics; applications we are most interested in determining how well the
382 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

TABLE III TABLE IV

ALGORITHM INFERENCE ABILITY ALGORITHM RECOVERY ABILITY

indicate that in this case, the estimation-exploration algorithm

caused the target robot to achieve the same performance as the
control algorithm, in which all evolution is performed directly
on the target robot, but using only 20 physical tests, compared
with 600. Thus, there was an order of magnitude improvement
in extracting a desired behavior from the target robot.
A similar control was devised for the damaged robot recovery
application. Again a genetic algorithm was used to evolve neural
controllers directly on the target robot, which in this case was
algorithm can approximate an indescribable structural change a damaged quadrupedal robot (one of the lower legs separates
to the robot or its environment and recover from it. completely from the body). In this case, even after 1550 evalu-
ations (evolving a population of 100 controllers for 30 genera-
D. Comparison of Performance tions), only half of the mean distance was achieved compared
with the mean distance traveled by the target robot after ten
Table III summarizes the inference ability of the algorithm
physical tests using the proposed algorithm. Although there is
for the grammar induction and gene network inference applica-
a slight overlap in the standard deviations, it is clear that the
tions. The first row reports that for 21 of the 30 target DFAs,
proposed algorithm performs better than simply evolving con-
the estimation-exploration algorithm was able to find a perfect
trollers directly on the target robot, using two orders of magni-
model using fewer target evaluations than the other two control
tude fewer target trials (1500 versus 10).
algorithms [see Fig. 7(b)]. The remaining nine targets showed
no statistical improvement using intelligent testing. Moreover,
for 14 of the 30 target DFAs, the proposed algorithm discovered VIII. CONCLUSION
a perfect model using significantly fewer model evaluations [see
Fig. 7(a)]. In the case of gene network inference, all three algo- We have introduced the estimation-exploration algorithm,
rithms performed 10 model evaluations, the standard GA a coevolutionary approach to system identification. We have
performed 100 target evaluations, the iterative algorithm with shown how our algorithm differs from other pure coevolu-
random testing performed 56 target evaluations and the estima- tionary systems because it does not so much introduce a new
tion-exploration algorithm performed 55 target evaluations. As advance in coevolutionary theory, but rather introduces a sys-
can be seen from Table III, the proposed algorithm consistently tematic, domain independent method for performing synthesis
produces more accurate models for large and more densely con- or analysis using coevolution.
nected gene networks (see Fig. 9). This indicates that the esti- Unlike other coevolutionary systems, the estimation-explo-
mation-exploration algorithm is scalable: intelligent testing be- ration algorithm evolves candidate models that are influenced
comes increasingly useful as the size of the target system in- not only by the competing population of tests, but also by the
creases and its observability decreases. target system. This removes two potential pathologies from the
Table IV summarizes the recovery ability of the proposed al- algorithm (cycling and overspecialization) which are often suf-
gorithm for the latter two robot applications. In the third appli- fered by other coevolutionary algorithms. The third pathology,
cation, the target robot had a different mass distribution and sig- disengagement, has been observed to occur in our algorithm, but
nificant sensor time lags, compared with its corresponding robot not for the four applications shown here. We have formulated
simulator. As a control, a neural controller was evolved directly two mechanisms that combat disengagement, and are currently
on this target robot using the genetic algorithm described in Sec- investigating additional mechanisms.
tion V. Ten independent runs of this control were performed. We have documented its application to four very different
After 600 evaluations in each run, the best controller was ex- problems, thus demonstrating its problem domain indepen-
tracted, and the mean distance achieved by the robot using these dence. The algorithm has several other benefits that have been
ten networks in turn was found to be 5.3 m. The intervals in highlighted in the various applications: coupled automation of
this table represent standard deviation, indicating the distribu- both models and intelligent tests; the active minimization of
tion of distances traveled by the best robot from each run, for target testing; increased usefulness as target systems increase
both algorithms. Both the means and the standard deviations in size and complexity; increased usefulness for target systems
BONGARD AND LIPSON: NONLINEAR SYSTEM IDENTIFICATION USING COEVOLUTION OF MODELS AND TESTS 383

with low observability; and automating the extraction of useful [7] J. Bongard and R. Pfeifer, “Repeated structure and dissociation of geno-
behavior from a hidden system. typic and phenotypic complexity in artificial ontogeny,” in Proc. Genetic
Evol. Comput. Conf., 2001, pp. 829–836.
Thus, our algorithm performs well at inference: it can au- [8] J. C. Bongard and H. Lipson, “Automated robot function recovery after
tomatically construct a hidden system that explains the gath- unanticipated failure or environmental change using a minimum of hard-
ered data. Active learning is also concerned with the intelli- ware trials,” in Proc. NASA/DoD Conf. Evolvable Hardware, Seattle,
WA, 2004, pp. 169–176.
gent selection of tests [16], but in that domain a classifier is [9] J. C. Bongard and H. Lipson, “Automating genetic network inference
used to describe input/output relations; no explicit model of the with minimal physical experimentation using coevolution,” in Proc. Ge-
target system is generated. This is an important distinction to netic Evol. Comput. Conf., Seattle, WA, 2004, pp. 333–345.
[10] D. W. Bradley and A. M. Tyrrell, “Immunotronics: Novel finite-state-
make, because by building explicit models to explain the data machine architectures with built-in self-test using self-nonself differen-
rather than learning an approximation of the input/output rela- tiation,” IEEE Trans. Evol. Comput., vol. 6, no. 3, pp. 227–238, Jun.
tionships, we obtain several benefits: 2002.
[11] S. Brnyjolfsson and A. Arnstrom, “Error detection and recovery in flex-
1) we can use the resultant model to learn more about the ible assembly systems,” Int. J Adv Manuf. Technol., vol. 5, pp. 112–125,
hidden system (we can use it to test hypotheses about the 1990.
[12] A. Bucci, J. B. Pollack, and E. D. De Jong, “Automated extraction of
actual system); problem structure,” in Proc. Genetic Evol. Comput. Conf., 2004, pp.
2) we can uncover possible causalities rather than just cor- 501–512.
relations (“gene regulates gene ,” rather than “the ex- [13] T. Chen, H. L. He, and G. M. Church, “Modeling gene expression with
differential equations,” in Proc. Pacific Symp. Biocomput., vol. 4, 1999,
pressions of genes and are correlated”); pp. 29–40.
3) we can obtain a comprehensible description of some [14] O. Cicchello and S. C. Kremer, “Inducing grammars from sparse data
change in the hidden system (such as what damage has sets: A survey of algorithms and results,” J. Mach. Learn. Res., vol. 4,
pp. 603–632, 2003.
incurred, and where); [15] D. Cliff and G. F. Miller, “Tracking the red queen: Measurements of
4) we can generate some useful behavior for it (continuation adaptive progress in coevolutionary simulations,” in Proc. Eur. Conf.
of function by circumventing the damage). Artif. Life, 1995, pp. 200–218.
[16] D. Cohn, L. Atlas, and R. Ladner, “Improving generalization with active
The third and fourth property allow us to recover function learning,” Mach. Learn., vol. 15, no. 2, pp. 201–221, 1994.
from the hidden system using a minimum number of tests. [17] H. de Jong, “Modeling and simulation of genetic regulatory systems: A
Second, unlike active learning, tests can be constructed, al- literature review,” J. Comput. Biol., vol. 9, no. 1, pp. 69–105, 2002.
[18] P. D’haeseleer, S. Liang, and R. Somogyi, “From co-expression
lowing for the exploration of a possibly infinite space of tests, clustering to reverse engineering,” Bioinformatics, vol. 16, no. 8, pp.
rather than selected from some fixed space of tests with equal 707–726, 2000.
dimensionality. [19] E. A. DiPaolo, “Homeostatic adaptation to inversion of the visual field
and other sensorimotor disruptions,” in From Animals to Animats 6, J.
Furthermore, using evolution allows us to synthesize models A. Meyer, A. Berthoz, D. Floreano, H. L. Roitblat, and S. W. Wilson,
and tests for systems about which little is known. This is impor- Eds. Cambridge, MA: MIT Press, 2000, pp. 440–449.
tant in nonlinear black-box system identification tasks for which [20] E. Z. Evans and S. G. Lee, “Automatic generation of error recovery
knowledge through learned activity,” in Proc. IEEE Int. Conf. Robotics
not even the underlying topology of the hidden system is known Autom., 1994, pp. 2915–2920.
[56]. Artificial evolution could be used to build explicit models [21] R. Featherstone and D. E. Orin, “Robot dynamics: Equations and algo-
directly from observed data in such instances, a challenge that rithms,” in Proc. IEEE Int. Conf. Robotics Autom., 2000, pp. 826–834.
[22] S. G. Ficici and J. B. Pollack, “Pareto optimality in coevolutionary
has not yet been addressed in the system identification literature, learning,” in Proc. 6th Eur. Conf. Adv. Artif. Life:, 2001, pp. 316–327.
to the best of the authors’ knowledge. For example the algorithm [23] D. Floreano and F. Mondada, “Evolutionary neurocontrollers for au-
could be used in a remote robotics application to synthesize a tonomous mobile robots,” Neural Netw., vol. 11, pp. 1461–1478, 1998.
[24] D. Floreano and J. Urzelai, “Neural morphogenesis, synaptic plasticity,
model of a novel environment based only on the robot’s sensor and evolution,” Theory Biosci., vol. 120, pp. 225–240, 2001.
data: this is an attractive avenue for future study. In future work, [25] L. J. Fogel, A. J. Owens, and M. J. Walsh, Artificial Intelligence Through
we also plan to apply the algorithm to a range of actual physical Simulated Evolution. New York: Wiley, 1966.
[26] P. Funes and J. Pollack, “Computer evolution of buildable objects,” in
systems, including physical robots and real biological networks. Evolutionary Design by Computer, P. Bentley, Ed. San Francisco, CA:
Morgan Kauffman, 1999, pp. 387–403.
[27] M. Gevers, “A decade of progress in iterative control design: From
REFERENCES
theory to practice,” J. Process Control, vol. 12, no. 4, pp. 519–531, May
[1] M. G. Abu-Hamdan and A. S. El-Gizawy, “Computer aided monitoring 2002.
system for flexible assembly operations,” Comput. Ind., vol. 34, pp. [28] J. J. Grefenstette and C. L. Ramsey, “An approach to anytime learning,”
1–10, 1997. in Proc. 9th Int. Workshop Mach. Learn., 1992, pp. 189–195.
[2] A. Adamatzky, M. Komosinski, and S. Ulatowski, “Software review: [29] W. D. Hillis, “Coevolving parasites improve simulated evolution as an
Framsticks,” Kybernetes: Int. J. Syst. Cybern., vol. 29, pp. 1344–1351, optimization procedure,” Physica D, vol. 42, pp. 228–234, 1990.
2000. [30] G. S. Hornby and J. B. Pollack, “Creating high-level components with
[3] A. Arkin, P. Shen, and J. Ross, “A test case for correlation metric con- a generative representation for body-brain evolution,” Artif. Life, vol. 8,
struction of a reaction pathway from measurements,” Science, vol. 277, no. 3, pp. 223–246, 2002.
pp. 1275–1279, 1997. [31] H. Iba and A. Mimura, “Inference of a gene regulatory network by means
[4] C. M. Baydar and K. Saitou, “Off-line error prediction, diagnosis and of interactive evolutionary computing,” Inf. Sci., vol. 145, pp. 225–236,
recovery using virtual assembly systems,” in Proc. IEEE Int. Conf. 2002.
Robotics Autom., 2001, pp. 818–823. [32] N. Jakobi, “Evolutionary robotics and the radical envelope of noise hy-
[5] F. Bergadano and D. Gunetti, Inductive Logic Programming: From Ma- pothesis,” Adaptive Behav., vol. 6, no. 1, pp. 131–174, 1997.
chine Learning to Software Engineering. Cambridge, MA: MIT Press, [33] E. D. De Jong and J. B. Pollack, “Ideal evaluation from coevolution,”
1995. Evol. Comput., vol. 12, no. 2, pp. 159–192, 2004.
[6] J. Bongard and H. Lipson, “‘Managed challenge’ alleviates disengage- [34] H. Juillé, “Methods for statistical inference: Extending the evolutionary
ment in coevolutionary system identification,” in Proc. Genetic Evol. computation paradigm,” Ph.D. dissertation, Brandeis Univ., Waltham,
Comput. Conf., 2005, pp. 531–538. MA, May 1999.
384 IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 9, NO. 4, AUGUST 2005

[35] J. F. Kao, “Optimal recovery strategies for manufacturing systems,” Eur. [57] S. Srinivas, “Error recovery in robot systems,” Ph.D. dissertation, Calif.
J. Oper. Res., vol. 80, pp. 252–263, 1995. Inst. Technol., Pasadena, CA, 1977.
[36] D. Keymeulen, M. Iwata, Y. Kuniyoshi, and T. Higuchi, “Online evolu- [58] K. O. Stanley and R. Miikkulainen, “Competitive coevolution through
tion for a self-adapting robotics navigation system using evolvable hard- evolutionary complexification,” J. Artif. Intell. Res., vol. 21, pp. 63–100,
ware,” Artif. Life, vol. 4, pp. 359–393, 1998. 2004.
[37] D. Keymeulen, A. Stoica, and R. Zebulum, “Fault-tolerant evolvable [59] A. Thompson, “Artificial evolution in the physical world,” in Evolu-
hardware using field programmable transistor arrays,” IEEE Trans. tionary Robotics: From Intelligent Robots to Artificial Life (ER’97), T.
Reliability (Special Issue on Fault-Tolerant VLSI Syst.), vol. 49, pp. Gomi, Ed. Ottawa, ON, Canada: AAI Books, 1997, pp. 101–125.
305–316, Sep. 2000. [60] S. Tokura, A. Ishiguro, H. Kawai, and P. Eggenberger, “The effect of
[38] S. Kikuchi, D. Tominaga, M. Arita, and M. Tomita, “Pathway finding neuromodulations on the adaptability of evolved neurocontrollers,” in
from given time-courses using genetic algorithm,” Genome Informatics, Proc. 6th Eur. Conf. Artif. Life, J. Kelemen and P. Sosik, Eds., 2001, pp.
vol. 12, pp. 304–305, 2001. 292–295.
[39] H. Kitano, Foundations of Systems Biology. Cambridge, MA: MIT [61] M. L. Visinsky, J. R. Cavallaro, and I. D. Walker, “Expert system frame-
Press, 2001. work for fault detection and fault tolerance in robotics,” Comput. Elect.
[40] J. Koza, Genetic Programming: On the Programming of Computers by Eng., vol. 20, pp. 421–435, 1994.
Means of Natural Selection. Cambridge, MA: MIT Press, 1992. [62] R. A. Watson and J. B. Pollack et al., “Coevolutionary dynamics in a
[41] S. Kumar and P. Bentley, On Growth, Form and Computers. New York: minimal substrate,” in Proc. Genetic Evolutionary Comput. Conf., L.
Academic, 2003. Spector and E. D. Goodman, Eds., 2001, pp. 702–709.
[42] NASA Jet Propulsion Laboratory. (2004) Mars exploration rovers. [On- [63] V. Zykov, J. Bongard, and H. Lipson, “Coevolutionary variance guides
line]. Available: http: //[Link]/mer2004/ physical experimentation in evolutionary system identification,” in Proc.
[43] H. Lipson and J. B. Pollack, “Automatic design and manufacture of ar- NASA/DoD Conf. Evolvable Hardware, 2005, pp. 213–220.
tificial lifeforms,” Nature, vol. 406, pp. 974–978, 2000. [64] S. G. Ficici, “Solution concepts in coevolutionary algorithms,” Ph.D.
[44] L. Ljung, System Identification: Theory for the User. Englewood dissertation, Comput. Sci. Dept., Brandeis Univ., Waltham, MA, 2004.
Cliffs, NJ: Prentice-Hall, 1999. Tech. Rep. CS-03-243.
[45] S. M. Lucas and T. J. Reynolds, “Learning DFA: Evolution versus evi-
dence driven state merging,” in Proc. Congr. Evol. Comput., 2003, pp.
351–358.
[46] S. Luke, S. Hamahashi, and H. Kitano, “‘Genetic’ programming,” in
Josh C. Bongard received the [Link]. degree (Honors)
Proc. Genetic Evol. Comput. Conf., 1999, pp. 1098–1105.
in computer science from McMaster University,
[47] S. H. Mahdavi and P. J. Bentley, “An evolutionary approach to damage
Hamilton, ON, Canada, in 1997, the [Link]. degree in
recovery of robot motion with muscles,” in Proc. 7th Eur. Conf. Artif.
evolutionary and adaptive systems from the School
Life, 2003, pp. 248–255.
of Cognitive and Computing Sciences, University
[48] E. Mjolsness, D. H. Sharp, and J. Reinitz, “A connectionist model of
of Sussex, Brighton, U.K. in 1999, and the Ph.D.
development,” J. Theor. Biol., vol. 152, pp. 429–454, 1991.
degree from the Artificial Intelligence Laboratory,
[49] B. Olsson, “Coevolutionary search in asymmetric spaces,” Inf. Sci., vol.
University of Sussex, for research in the field of
133, pp. 103–125, 2001.
evolutionary robotics.
[50] R. Parekh and V. Honavar, “Automata induction, grammar inference,
He is currently working as a Postdoctoral Re-
and language acquisition,” in The Handbook of Natural Language Pro-
searcher in the Sibley School of Mechanical and
cessing, R. Dale, H. Moisl, and H. Somers, Eds. New York: Marcel
Aerospace Engineering, Cornell University, Ithaca, NY.
Dekker, 2000, pp. 727–764.
[51] J. B. Pollack, H. Lipson, S. Ficici, P. Funes, G. Hornby, and R. Watson,
“Evolutionary techniques in physical robotics,” in Evolvable Systems:
From Biology to Hardware, J. Miller, Ed. Berlin, Germany: Springer-
Verlag, 2000, pp. 175–186. Hod Lipson received the [Link]. and Ph.D. degrees in
[52] M. A. Potter and K. A. De Jong, “Cooperative coevolution: An architec- mechanical engineering and computer-aided design
ture for evolving coadapted subcomponents,” Evol. Comput., vol. 8, no. from the Technion—Israel Institute of Technology,
1, pp. 1–29, 2000. Haifa.
[53] T. Reil and P. Husbands, “Evolution of central pattern generators for He joined the faculty of the Departments of
bipedal walking in a real-time physics environment,” IEEE Tran. Evol. Mechanical and Aerospace Engineering and Com-
Comput., vol. 6, no. 2, pp. 159–168, Apr. 2002. puting and Information Science, Cornell University,
[54] C. D. Rosin and R. K. Belew, “New methods for competitive coevolu- Ithaca, NY, in 2001. Prior to this appointment, he
tion,” Evol. Comput., vol. 5, no. 1, pp. 1–29, 1997. was a Postdoctoral Researcher in the Computer
[55] K. Sims, “Evolving 3D morphology and behavior by competition,” Artif. Science Department, Brandeis University and a
Life IV, pp. 28–39, 1994. Lecturer in the Mechanical Engineering Department,
[56] J. Sjöberg, Q. Zhang, L. Ljung, A. Benveniste, B. Delyon, P.-Y. Glo- Massachusetts Institute of Technology (MIT), Cambridge. Before joining
rennec, H. Hjalmarsson, and A. Juditsky, “Nonlinear black-box mod- academia, he spent several years as a research engineer in the mechanical,
eling in system identification: A unified overview,” Automatica, vol. 31, electronic, and software industries. His interests are in the area of design
no. 12, pp. 1691–1724, 1995. automation and evolutionary robotics.