\chapter{Genetic and Evolutionary Algorithms}
\begin{flushright}
\parbox{5in}{The word ``algorithm''\index{algorithm} itself is quite interesting; at first glance it may look as though someone intended to write ``logarithm'' but jumbled up the first four letters. The word did not appear in {\it Webster's New World Dictionary\/}\index{Webster's New World Dictionary} as late as 1957; we find only the older form ``algorism'' with its ancient meaning, i.e., the process of doing arithmetic using Arabic numerals.''$\ldots$``Finally, historians of mathematics found the true origin of the word algorism: it comes from the name of a famous Persian textbook author, Abu Ja`far Mohammed ibn M\^{u}s\^{a} al-Khow\^{a}rizm\^{\i} (c.~825)---literally, ``Father of Ja`far, Mohammed, son of Moses, native of Khow\^{a}rizm.'' ''$\ldots$``Al-Khow\^{a}rizm\^{\i} wrote the celebrated book {\it Kitab al jabr w'al-muqabala} (``Rules of restoration and reduction''); another word, ``algebra'', stems from the title of his book, although the book wasn't really very algebraic.''$\ldots$``The modern meaning for algorithm is quite similar to that of {\it recipe\/}, {\it process\/}, {\it method\/}, {\it technique\/}, {\it procedure\/}, {\it routine\/}, except that the word ``algorithm'' connotes something just a little different.~\cite[pages~1--3]{knuth73}}
\end{flushright}

\bigskip

Genetic algorithms,\index{genetic algorithm (GA)}\index{GA|see{genetic algorithm}} introduced by Holland~\cite{holland92} in the 1960s, have been originally proposed to model adaptive processes, both natural and man-made. In this chapter, we discuss applications of genetic algorithms to solving optimization problems. Genetic algorithms were inspired by biological genetics\index{biological genetics} and evolutionary processes.\index{evolutionary processes} This is the reason for biological terms being used in the context of genetic algorithms. One of the variants of genetic algorithms, called the evolutionary algorithm\index{evolutionary algorithm (EA)}\index{EA|see{evolutionary algorithm}} is also presented.

\section{Genetics as an Inspiration for an Optimization Approach}
The {\em gene\/} is the basic unit of inheritance.\index{gene!basic unit of inheritance} Genes are pieces of DNA,\index{genes!pieces of DNA} which is the material inside of a nucleus of a cell that carries the genetic instructions for making living organisms. Many traits are determined by genes inherited from the parents of the individual. Consider gene ${\cal A}$ that determines a certain trait, such as seed coat texture---round or wrinkled. Thus gene ${\cal A}$ has two alternative forms. These different forms of the same gene are called {\em alleles\/}, or {\em allelic genes}.\index{alleles!different forms of the same gene}\index{allelic genes|see{alleles}} The allelic forms of each gene are symbolized by the uppercase and lowercase of the same letter. In our example, gene ${\cal A}$ has two alleles, allele $A$ and allele $a$. One of the alleles controls the round seed coat texture while, the other allele controls the wrinkled texture. Genes are organized into {\em chromosomes},\index{chromosomes}\index{genes!organized into chromosomes} which are small packages of genes in the nucleus of a cell. Different kinds of organisms have different numbers of chromosomes. Humans have 23 pairs of chromosomes including two sex chromosomes. In mating each parent contributes one chromosome to each pair of the chromosomes of an offspring, so the offspring gets half of his or her chromosomes from his or her mother and half from his or her father. The assemblage of genes, encrypted as a string of DNA molecules, is called the {\em genome}.\index{genome} The human genome consists of about 35,000 genes. The particular gene constitution is called the {\em genotype}.\index{genotype!gene constitution} The expression of the genotype is the {\em phenotype\/}, which is the way we look and behave.\index{phenotype!expression of genotype}

Genetic algorithms are a result of an effort to model adaptation phenomena in natural and artificial systems. The term {\em adaptation\/} can be viewed as the composition of two words: {\em ad}~+~{\em aptare\/} meaning ``to fit to''~\cite[page~1]{holland92}. ``Just what are adaptation's salient features? We can see at once that adaptation, whatever its context, involves a progressive modification of some structure or structures.\index{adaptation!progressive modification of structures} These structures constitute the grist of the adaptive process, being largely determined by the field of study. Careful observation of successive structural modifications generally reveals a basic set of structural modifiers or operators; repeated action of these operators yields the observed modification sequences.''~\cite[page~3]{holland92}

In this chapter, we describe two genetic algorithms and apply them to solving optimization problems. An optimization problem of a real-valued function, $f:\Omega \to \mb{\Bb R}$, requires finding an element $\b x\in \Omega$ such that $f$ is maximized or minimized. The function $f$ is called the {\em fitness function\/}.\index{fitness function} Suppose that we wish to maximize $f$ over the parameter space $\Omega$.\index{parameter space} The parameter space, $\Omega$, of an underlying optimization problem is also called the {\em phenotype space\/}.\index{phenotype space}\index{parameter space!as phenotype space} A classical genetic algorithm works with binary representations of the points from the parameter space, rather than directly with the elements of $\Omega$. We refer to the classical genetic algorithms as the canonical algorithms.\index{canonical GA!classical GA} The strings of binary digits (1s and 0s) corresponding to these points are called chromosomes.\index{chromosomes!in canonical GA strings of binary digits} The space of chromosomes is called the {\em genotype space\/}.\index{genotype space!space of chromosomes} In general, the chromosomes in  the genotype space are abstract mathematical objects on which the so called genetic operators work. After a stopping criterion for the algorithm is satisfied, a decoding function transforms a solution chromosome from the genotype space into the phenotype space.  Genetic algorithms are derivative-free optimization algorithms\index{genetic algorithms!derivative-free optimization algorithms} that use only the fitness function. They do not use the fitness function derivatives. The higher the fitness of a chromosome, the better the candidate solution encoded by the chromosome. Genetic algorithms work with populations of chromosomes\index{genetic algorithms!work with populations of chromosomes} rather than individual chromosomes. A genetic algorithm is an iterative process that at each iteration uses the information obtained from testing a current population of chromosomes to direct its search into promising areas of the search space.\index{search space} Thus, generating a suitable population of chromosomes at each iteration is the key to the algorithm performance. In the simplest form, creating a new population of chromosomes to be tested involves three genetic operators: selection, crossover, and mutation. These genetic operators\index{genetic operators!work on genotypes} work on genotypes that encode phenotypes.\index{genotypes!encode phenotypes}

The role of the selection operator\index{genetic operators!selection operator}\index{selection operator} is to pick a population of chromosomes for further processing, called reproduction. Chromosomes are selected into a mating pool\index{mating pool} for reproduction with a probability proportional to their fitness.

The crossover operator\index{genetic operators!crossover operator}\index{crossover operator} randomly exchanges substrings between two parent chromosomes to create two offspring. The simplest crossover operator is the one-point crossover, where a crossing site is randomly picked and the corresponding substrings are exchanged, as illustrated in Figure~\ref{fig:oneptcr}. For example, the two parent chromosomes
\[
1010101\quad\mb{and}\quad 1111010,
\]
could have the crossing site after the fourth bit in each to produce the two offspring
\[
1010010\quad\mb{and}\quad 1111101.
\]
\begin{figure}
\centering
\includegraphics[scale=0.3]{oneptcr.eps}
%\vspace{3.5in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/oneptcr.bmp}
%\hskip1in\special{bmp:c:/mybook/oneptcr.bmp x=4.67in y=3.5in}
\caption{The one-point crossover.}
\label{fig:oneptcr}
\end{figure}

The mutation operator\index{genetic operators!mutation operator}\index{mutation operator} randomly flips an arbitrary bit in a chromosome. For example, the chromosome 010101 might undergo a mutation in its third position to yield a new chromosome 011101.

We now provide a description of the canonical genetic algorithm using the above components.
\begin{algorithm}\rm
\begin{center}
{\bf Canonical Genetic (GA)}
\end{center}
For a given optimization problem characterized by a fitness function $f$:\\
{\bf Step 1.} Choose a method of encoding candidate solutions into chromosomes, that is, decide upon a map between the phenotype and genotype spaces. (A commonly used encoding method is to concatenate the binary strings representing each real parameter value $x_i$.)\\
{\bf Step 2.} Start with a randomly generated population of chromosomes.\\
{\bf Step 3.} Evaluate the fitness of each chromosome in the population.\\
{\bf Step 4.} Replace the current population with a new population by applying selection, crossover, and mutation operations on the current population.\\
{\bf Step 5.} If a stopping criterion is satisfied, stop and choose the fittest chromosome from the population as a solution, else go to {\bf Step 3.}\end{algorithm}\markex
Another type of a genetic algorithm that we analyze in this chapter attempts to avoid the need for a decoding function by using the chromosome representation as close as possible to their representation in the phenotype space.

\section{Implementing a Canonical Genetic Algorithm}
We now discuss implementation issues of the canonical genetic algorithm~(GA)\index{canonical genetic algorithm (GA)} when applied to solving optimization problems. We consider optimization problems of the form:
\begin{equation}
\mb{maximize}\quad g(\b x),\quad \b x\in \Omega \subset\mb{\Bb R}^n.
\label{eq:maxproopt}
\end{equation}
A maximizing point or maximizer is denoted $\b x^*$. We note that there is no loss of generality by considering only maximization problems, for if we are asked to minimize $g$, we can transform the minimization problem into the maximization problem by means of a simple transformation,
\begin{equation}
\min_{\b x} g(\b x)=-\max_{\b x}(-g(\b x)).
\label{eq:mintomat}
\end{equation}
Also note that if $\b x^*$ is a minimizer of $g(\b x)$, then $\b x^*$ is a maximizer of $-g(\b x)$, and vice versa. Once an objective function of the maximization problem is given, the next step is to transform the objective function into a fitness function.\index{objective function!transforming into fitness function}\index{fitness function} For example, we can use the transformation,
\begin{equation}
f(\b x)=g(\b x)-g_{\min},
\label{eq:conobjifi}
\end{equation}
where $f(\b x)$ denotes a fitness function and $g_{\min}$ is the minimum value of $g(\b x)$ for the current population of chromosomes. Note that the constructed fitness function is indeed non-negative on the chromosomes of the current population.

We now discuss the implementation of {\bf Step 1} of the canonical GA presented in the previous section. We first analyze the case when the fitness function, $f$, is a function of only one parameter $x$. Suppose that
\[
x\in \left[x_{\min}, x_{\max}\right]=\Omega.
\]
A common encoding method, that transforms points from the phenotype space into a genotype space, is the binary encoding.\index{binary encoding!phenotype space into a genotype space} In our example the phenotype space is the interval $\left[x_{\min} ,x_{\max}\right]$. If we decide to use $l$ bits to represent a point $x$, then this corresponds to mapping the interval $\left[x_{\min},x_{\max}\right]$ into the discrete genotype space\index{genotype space!discrete} $\left[0,2^l-1\right]$ consisting of $2^l$ binary strings. Thus, the number of bits $l$ determines the size of the genotype space and hence the resolution, or accuracy, of the subsequent search. We can control the resolution by appropriately choosing the parameter $l$. Let $r$ denote the desired resolution. Then, the parameter $l$ that results in the required resolution can be found by solving the inequality
\begin{equation}
r\ge\frac{x_{\max}-x_{\min}}{2^l-1}.
\label{eq:desresga}
\end{equation}

If the fitness function depends on many parameters, then we can encode the elements of the phenotype space by concatenating the binary strings representing each component of $\b x=\left[\begin{array}{cccc}
x_1 & x_2 & \cdots & x_n
\end{array}\right]^{\top}$. For example, let $f=f(x_1,x_2)$, where $x_1\in [x_{1,\min},x_{1,\max}]$ and $x_2\in [x_{2,\min},x_{2,\max}]$. If we use 10 bits to encode each variable, then the genotype space will consist of 20-bit binary strings and the size of the genotype space is $2^{10}\times 2^{10}=1,048,576$.

After initializing randomly a population of chromosomes, we evaluate the fitness $f(\b x)$ of each chromosome $\b x$ in the population. This requires decoding from the genotype space into the phenotype space.\index{decoding!from genotype space into phenotype space} Suppose, as before, that the fitness function depends on two parameters $x_1$ and $x_2$ and that we used the standard encoding method described above. We now have to convert a given binary string into two real numbers $x_1\in [x_{1,\min},x_{1,\max}]$ and $x_2\in [x_{2,\min},x_{2,\max}]$. We can perform the required decoding process in three steps. First note that the first $l$ bits of the chromosome $\b x$ correspond to $x_1$ while the remaining bits correspond to $x_2$. We first convert the strings of $l$ bits into integers. Thus, for example, if $l=4$ and the chromosome is 10010111, which is an 8-bit string, then the first 4 bits yield
\[
x_{1,int}=1\times 2^3+0\times 2^2+ 0\times 2^1 +1\times 2^0=9,
\]
while the remaining 4 bits yield
\[
x_{2,int}=0\times 2^3+1\times 2^2+ 1\times 2^1 +1\times 2^0=7.
\]
Then, we convert $x_{1,int}$ from the range $[0, 2^l-1]$ into the range $\left[0,\left(x_{1,\max}-x_{1,\min}\right)\right]$ by multiplying $x_{1,int}$ by
\begin{equation}
\frac{x_{1,\max}-x_{1,\min}}{2^l-1}.
\label{eq:movranxo}
\end{equation}
To obtain the corresponding value of $x_1$, we add $x_{1,\min}$ to~(\ref{eq:movranxo}) to get
\begin{equation}
x_1=x_{1,\min}+x_{1,int}\frac{x_{1,\max}-x_{1,\min}}{2^l-1}.
\label{eq:voxoings}
\end{equation}
We obtain $x_2$ proceeding in an analogous manner. In general, if $\b x$ has $n$ components, the mapping from the genotype space into the phenotype space\index{genotype space!into phenotype space}\index{phenotype space!from genotype space} has the form
\begin{equation}
\fbox{$\displaystyle x_i=x_{i,\min}+x_{i,int}\frac{x_{i,\max}-x_{i,\min}}{2^l-1},\quad i=1,2,\ldots,n $}
\label{eq:genxoings}
\end{equation}

The next step of the GA algorithm is the application of the selection operator to pick chromosomes from the current population for reproduction to create the new population. In the canonical GA, a pair of chromosomes is selected from the current population with the probability of selection being proportional to the fitness values of the chromosomes. Selection is done with replacement, which means that the same chromosome can be selected more than once for reproduction. A common method of selecting parent chromosomes\index{selecting parent chromosomes!roulette-wheel method} is the so called {\em roulette-wheel method}.\index{roulette-wheel method!selecting parent chromosomes} In this method each chromosome is given a slice of a wheel. The area of the wheel slice is proportional to the chromosome's fitness. The wheel is spun, and the chromosome is selected on whose wheel slice the ball comes to rest. The roulette-wheel method can be implemented in the following steps:
\begin{enumerate}
\item Calculate the fitness of each chromosome in the current population;
\item Calculate the total fitness of the current population by summing the fitnesses of each chromosome in the current population;
\item Divide the fitness of each chromosome by the total fitness;
\item \label{item:cumfit} Calculate the cumulative fitness,\index{cumulative fitness} labeled for our purposes as \verb!cum_fitness!, of the current population using the scaled fitness of each chromosome found in the previous step. You can use MATLAB's command \verb!cumsum! for this purpose.
\item Return the first population member whose corresponding cumulative fitness is greater than the random number from the interval $(0,1)$. The following MATLAB command can be used here: \verb!find(cum_fitness-rand>0)!, where \verb!cum_fitness! is the vector of the scaled cumulative fitnesses found in step~\ref{item:cumfit}.
\item Repeat the above step to get the second parent chromosome.
\end{enumerate}

Using the parent chromosomes picked by the selection operator,\index{selection operator} we form two offspring\index{offspring} by crossing over, with probability $p_c$, at a randomly chosen crossing site.\index{crossing site} The uniform probability $p_c$ is called the crossover probability\index{crossover probability} or crossover rate.\index{crossover rate} If no crossover takes place, two offspring chromosomes are formed that are exact copies of their parent chromosomes. One can also use a multi-point crossover operator.\index{crossover operator!multi-point}\index{multi-point!crossover operator} The two-point crossover operation is illustrated in Figure~\ref{fig:twoptcr}.\index{crossover operator!two-point}\index{two-point!crossover operator}
\begin{figure}
\centering
\includegraphics[scale=0.3]{twoptcr.eps}
%\vspace{3.5in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/twoptcr.bmp}
%\hskip.7in\special{bmp:c:/mybook/twoptcr.bmp x=4.67in y=3.5in}
\caption{The two-point crossover.}
\label{fig:twoptcr}
\end{figure}

The offspring chromosomes\index{offspring chromosomes} obtained from the crossover, are then acted upon by the mutation operator. Mutation is performed at each bit of the offspring chromosomes with probability $p_m$. The probability $p_m$ is called the mutation probability\index{mutation probability} or mutation rate.\index{mutation rate} The mutation operator randomly flips some of the bits in a chromosome. For example, the chromosome 101010 might be mutated in its third bit to yield 100010. The mutation probability is usually very small, for example $p_m=0.01$.

Each iteration, also called a {\em generation\/},\index{GA generation!iteration of GA} of the GA replaces the current population with a new one. A typical number of iterations of a GA algorithm is from 50 to 500 or more. The entire set of generations is called a {\em run\/}.\index{GA run!entire set of generations}
\begin{example}\rm
We apply the canonical GA described above to maximize the function,
\begin{equation}
f(x)=-15\sin^2(2x)-(x-2)^2+160,\quad |x|\le 10.
\label{eq:manminfu}
\end{equation}
This function has many local maximizers and minimizers as can be seen from its plot shown in Figure~\ref{fig:manminfu}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{manminfu.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/mybook/manminfu.bmp}
%\hskip1.2in\special{bmp:c:/mybook/manminfu.bmp x=3.97in y=3.5in}
\caption{A plot of the function given by~(\ref{eq:manminfu}).}
\label{fig:manminfu}
\end{figure}
The canonical GA was applied to maximize the function given by~(\ref{eq:manminfu}) with the chromosomes composed of 12 bits, $p_c=0.75$, and $p_m=0.01$. In Figure~\ref{fig:onegenr} a summary of a typical run is shown. In this run the best candidate solution was $x=1.565324$ whose corresponding objective value is $f(x)=159.809259$. The best solution that was observed in other runs resulted in $f(1.575092)=159.818346$.
\begin{figure}[h]
\centering
\includegraphics[scale=0.55]{onegenr.eps}
%\vspace{4in}
%\centerbmp{4.54in}{4in}{c:/mybook/onegenr.bmp}
%\hskip.75in\special{bmp:c:/mybook/onegenr.bmp x=4.54in y=4in}
\caption{Summary of a sample run of the canonical GA when searching for a maximizer of the function given by~(\ref{eq:manminfu}).}
\label{fig:onegenr}
\end{figure}
\end{example}\markex

In the next section, we attempt to obtain answers to some whys and hows of genetic algorithms.

\section{Analysis of the Canonical Genetic Algorithm}

In the analysis of genetic algorithms, the notion of a {\em schema\/}\index{schema!a set of bit strings} is useful. Plural of ``schema'' is {\em schemas\/}\index{schemas!plural of ``schema''} or {\em schemata\/}.\index{schemata!plural of ``schema''} A schema is a set of bit strings. It can be described using a template made up of ones, zeros, and asterisks that represent ``don't cares''. Thus, the asterisk~$\ast$ matches either 0 or 1. For example, the schema $H=\ast 0 \ast $ contains the binary strings $0 0 0$, $0 0 1$, $1 0 0$, and $1 0 1$. Note that a bit string of length $l$ belongs to $2^l$ different schemas. For example, the string $1 0$ of length $l=2$ belongs to the schemas $\ast \ast$, $\ast 0$, $1 \ast$, and $1 0$. Similarly, when $l=3$, then any bit string of length $l=3$ belongs to $2^3=8$ schemas. Elements of a set of $N$ bit strings of length $l$ are members of between $2^l$ and $N\times 2^l$ schemas. If all the strings are identical, then they belong to $2^l$ schemas, otherwise there may be up to $N\times 2^l$ schemas to which the $N$ strings of length $l$ can belong.

The parameters that characterize a given schema are its {\em order\/} and {\em defining length\/}.\index{order!of schema}\index{defining length!of schema} The order of a schema $H$, denoted $o(H)$, is the number of its fixed bits.\index{order!number of schema's fixed bits} For example, the order of $H=\ast 0 1\ast 1$ is $o(H)=3$ because $H$ has three bits fixed. The defining length of the schema $H$, denoted $d(H)$, is the distance between its first and the last fixed bits,\index{defining length!distance between first and last fixed bits} which is the difference between the positions of the rightmost fixed symbol and the leftmost fixed symbol. The defining length of $H=\ast 0 \ast 1$ is $d(H)=4-2=2$. In further discussion, the term schema is used to denote a set of binary strings represented by a template, as well as the template itself.\index{schema!as template}

The capability of a chromosome to solve a given optimization problem, that is, the chromosome performance, is measured by its fitness.\index{chromosome fitness}\index{fitness!of chromosome}\index{chromosome performance!its fitness} Similarly, we can characterize the performance of a schema by its average fitness.\index{schema performance!its average fitness} We now derive an expression for the average fitness of a given schema. Let $\b x$ be a chromosome that matches the template $H$, that is, $\b x \in H$. Then, we say that ``$\b x$ is an instance of $H$.'' Let $H$ be a schema with at least one instant in the $k$-th generation. Denote by $e(H,k)$ the number of instances of $H$ in the $k$-th generation, and let $\bar{u}(H,k)$ be the average fitness of $H$ in the $k$-th generation, that is,
\begin{equation}
\bar{u}(H,k)=\frac{\sum_{\b x\in H}f(\b x)}{e(H,k)}.
\label{eq:avefitohs}
\end{equation}
When the genetic algorithm is explicitly evaluating the fitness of $N$ chromosomes from the current generation, it is implicitly evaluating the average fitness of up to $N\times 2^l$ schemas to which the $N$ chromosomes belong.

We now find ${\cal E}\left(H,k+1  \right)$, the expected number of instances of the schema $H$ in the $(k+1)$-st generation, given the $k$-th generation. We first analyze the effect of the selection operator\index{selection operator} only. Assume that the roulette-wheel method is used in selecting parent chromosomes to the mating pool for reproduction. The probability that the chromosome $\b x$ is selected to the mating pool, $M(k)$, is
\[
\frac{f(\b x)}{\sum_{i=1}^N f\left(\b x_i^{(k)}\right)},
\]
where $\b x_i^{(k)}$ refers to the $i$-th chromosome in the $k$-th generation and  $f(\b x)$ is the fitness of the chromosome $\b x$. Therefore, the expected number of chromosomes equal to $\b x$ that are selected to the mating pool is
\[
N\frac{f(\b x)}{\sum_{i=1}^N f\left(\b x_i^{(k)}\right)}=\frac{f(\b x)}{\bar{f}(k)},
\]
where
\[
\bar{f}(k)=\frac{1}{N} \sum_{i=1}^N f\left(\b x_i^{(k)}\right)
\]
is the average fitness of the $k$-th generation.\index{average fitness!of $k$-th generation} Let $m(k)$ be the number of chromosomes in $M(k)$ that are instances of $H$. The expected value, ${\cal M}(H,k)$, of the number of chromosomes in the mating pool that are instances of $H$ is
\begin{equation}
\frac{\sum_{\b x\in H}f(\b x)}{\bar{f}(k)} = \frac{\sum_{\b x\in H}f(\b x)}{\bar{f}(k)}\frac{e(H,k)}{e(H,k)}.
\label{eq:avefkpogk}
\end{equation}
Substituting~(\ref{eq:avefitohs}) into~(\ref{eq:avefkpogk}) yields an alternative expression for the expected value of the number of chromosomes in the mating pool that are instances of $H$,
\begin{equation}
\fbox{$\displaystyle {\cal M}(H,k)=\frac{\bar{u}\left( H, k \right)}{\bar{f}(k)} e\left(H, k \right) $}
\label{eq:avefkpofi}
\end{equation}
It follows from the above equation that a schema whose average fitness is greater than the population average, will be found in an increasing number of chromosomes.

We next quantify the effects of the one-point crossover\index{crossover!one-point} on a given schema. Given two chromosomes, where at least one of them is an instance of the schema $H$, we say that the schema $H$ survives crossover\index{schema!survives crossover} operation if at least one of the offspring is an instance of $H$. We will give a lower bound on the probability $S_c(H)$ that the schema $H$ will survive the one-point crossover. A schema survives if the crossing site occurs outside the defining length. Crossovers that take place within the defining length can destroy the schema. For example, let $H=\ast 0 0 \ast \ast$, and let the parent chromosomes be $p_1=10011$ and $p_2=11100$. Note that $p_1\in H$, while $p_2 \notin H$. For the crossing site after the third position, there are two offspring chromosomes: $s_1=10000$, which is an instance of $H$, and $s_2=11111$, which is not an instance of $H$. Thus, the schema $H$ survived this crossover. Suppose now that the crossing site is after the second position, which is within the defining length of the schema $H$. The offspring chromosomes are $s_1=10100$ and $s_2=11011$. None of the offspring chromosome is an instance of $H$, which means that $H$ was destroyed by the crossover. If both parent chromosomes are instances of $H$, then the offspring chromosomes are also instances of $H$ irrespective of the location of the crossing site. An upper bound on the probability $D_c(H)$ that the schema $H$ will be destroyed by the one-point crossover is
\begin{equation}
D_c(H)\le p_c\frac{d(H)}{l-1}.
\label{eq:descrba}
\end{equation}
Therefore, a lower bound on the probability $S_c(H)$ that $H$ survives is
\begin{equation}
\fbox{$\displaystyle S_c(H)= 1-D_c(H)\ge 1-p_c\frac{d(H)}{l-1} $}
\label{lobosur}
\end{equation}
Thus, the probability of survival under one-point crossover is higher for short schemas.\index{schema!probability of survival}

A schema $H$ survives under mutation of its element, instance, if all defined bits of the instance of $H$ remain unchanged. Let $S_m(H)$ denote the probability that the schema $H$ survives under mutation of an instance of $H$. For each defined bit, the probability that it will not be mutated is $1-p_m$. The probability that no defined bit of $H$ will be mutated is obtained by multiplying $1-p_m$ by itself $o(H)$ times. Therefore,
\begin{equation}
\fbox{$\displaystyle S_m(H)=\left(1-p_m \right)^{o(H)} $}
\label{eq:mutsurp}
\end{equation}

We combine~(\ref{eq:avefkpofi}), (\ref{lobosur}), and~(\ref{eq:mutsurp}) to obtain a lower bound on the expected number of instances of $H$ in the $(k+1)$-st generation, given the $k$-generation,
\begin{equation}
\fbox{$\displaystyle {\cal E}\left(H,k+1 \right)\ge\frac{\bar{u}\left( H, k \right)}{\bar{f}(k)} e\left(H, k \right)\left(1-p_c\frac{d(H)}{l-1} \right) \left(1-p_m \right)^{o(H)} $}
\label{eq:schtheor}
\end{equation}
The above is known as the {\em schema theorem\/} of Holland.\index{schema theorem of Holland}\index{Holland's schema theorem} It follows from our discussion that the canonical genetic algorithm implicitly processes, in parallel, schemas. This implicit computation is performed with no operation being explicitly performed on schemas. To illustrate the last statement, consider the string $1 0 1$. As the fitness of this string is evaluated, partial information is obtained about the expected fitness of all possible schemas to which the evaluated string belongs. In our example, computing the fitness of the string $1 0 1$ gives partial information about the fitness of strings in the schemas $\ast \ast \ast $, $\ast \ast 1$, $\ast 0 \ast $, $1 \ast  \ast$, $\ast 0 1$, $1\ast 1$,\ldots\,.

\section{Simple Evolutionary Algorithm}
\label{sec:simpleea}
Canonical genetic algorithms use binary representations\index{canonical GAs!use binary representations} of the elements of the parameter space. The reason behind this is the schema theory, discussed in the previous section, which allows one to analyze canonical genetic algorithms. When using a binary representation of the elements of the parameter space, canonical algorithms use encoding and decoding functions. A binary representation may have some disadvantages though. Indeed, let $h:\{0, 1\}^l \to \Omega$ denote a decoding function used in the canonical genetic algorithm implementation. Then the original problem of optimizing $f:\Omega \to \mb{\Bb R}$ is transformed into the problem of optimizing the combined objective function
\[
f(h): \{0, 1\}^l\to\mb{\Bb R}.
\]
The combined objective function to be optimized may be more complex than the original objective function. For example, it may introduce extra nonlinearities or add extra optimizers, which can make the search process for a global optimizer\index{global optimizer} more difficult.

In contrast to the canonical genetic algorithms, in the optimization procedure discussed next, the parameter space is no longer coded using the binary representation. The genetic operators work directly on real-valued vectors. A term used to describe an optimization algorithm using genetic-type operators working on the real-valued vectors is an {\em evolutionary algorithm\/}.\index{evolutionary algorithm (EA)}\index{EA|see{evolutionary algorithm (EA)}} In what follows, we present a very simple evolutionary algorithm for solving optimization problems of the form given by~(\ref{eq:maxproopt}).
\begin{algorithm}\rm
\begin{center}
{\bf Simple Evolutionary Algorithm (EA)}
\end{center}
For a given optimization problem of the form~(\ref{eq:maxproopt}) perform the following steps:\\
{\bf Step 1.} Initialize a population of candidate solutions (chromosomes) to the problem. Each chromosome is a real-valued vector of the same length as the vector $\b x$ of decision variables.\\
{\bf Step 2.} Evaluate the fitness of each candidate solution in the population.\\
{\bf Step 3.} Replace the current population with a new population by applying selection, crossover, and mutation operations on the current population. The roulette-wheel selection method can be used for selection. A possible implementation of the crossover operation is to use the average crossover operator. This operator takes two parents, say $\b x_j$ and $\b x_k$ and returns one offspring
\[
\b x=\frac{1}{2}\left( \b x_j+\b x_k\right).
\]
Another implementation of the crossover operator is a convex combination of two parents. Specifically, let $r\in (0,1)$. Then, the offspring is
\[
\b x=r\b x_j+(1-r)\b x_k .
\]
Note that if the feasible set is convex, then the above crossover operators produce chromosomes that are elements of the feasible set. Otherwise conditions should be imposed to ensure that candidate solutions are the elements of the feasible set.

As for the mutation operation, let $\b x$ be a chromosome to be mutated. Then, an implementation of the mutation operation could have the form,
\[
\b x^{\prime}=\b x+\b d,
\]
where $\b d$ is a randomly generated feasible vector so that $\b x^{\prime}$ is in the feasible set.\\
{\bf Step 4.} If a stopping criterion is satisfied, terminate the algorithm and choose the fittest chromosome from the population as a solution, else go to {\bf Step 2.}\end{algorithm}\markex


In addition, both in the canonical GA, as well as in the evolutionary algorithm, the elitist strategy,\index{elitist strategy!introduced by K. DeJong} introduced by K.~DeJong in 1975, can be incorporated. This strategy copies the best chromosomes from the old population into the next population to prevent losing the best solutions in the succeeding iterations. As observed by Mitchell~\cite{mitchell96}, among others, the elitist strategy may lead to the domination of a population by a ``super chromosome''. However, the elitist strategy improves genetic and evolutionary algorithms performance. Rudolph~\cite{rudolph94} showed in 1994 that a canonical GA will never converge to the global optimum unless the best solution in the population, either before or after selection, is copied into the next population.

\section{Evolutionary Fuzzy Logic Controllers}

In this section, we employ evolutionary algorithms to design tracking controllers.\index{evolutionary algorithms!to design tracking controllers} We describe three design techniques for evolutionary fuzzy logic controllers\index{evolutionary fuzzy logic controllers} for the step-lane-change maneuver of a ground vehicle.\index{step-lane-change maneuver!of a ground vehicle} The presentation of this section is based on~\cite{yonggon02}. We mention that the step-lane-change maneuver is an important aspect in the design of the Intelligent Vehicle/Highway System~(IVHS),\index{Intelligent Vehicle/Highway System (IVHS)} the highway of the future.


An essential element in combining an evolutionary algorithm~(EA)\index{evolutionary algorithm (EA)} and a fuzzy logic controller~(FLC)\index{fuzzy logic controller (FLC)} is the encoding method used to represent a fuzzy rule base as a chromosome on which the EA operates. Selecting an appropriate encoding procedure is a critical factor that affects the efficiency of the EA and the final solution. Another essential element is the fitness function. Each chromosome defines a fuzzy rule base. After
conducting a simulation of the structure depicted in Figure~\ref{fig:gfsys} and measuring the tracking error, the fitness of the chromosome representing the fuzzy rule base of the simulated FLC is calculated.
\begin{figure}
\centering
\includegraphics[scale=0.75]{Gfsys2.eps}
%\vspace{4.25cm}
%\hskip1.15in\special{eps:c:/mybook/Gfsys2.eps y=4.25cm}
\caption{A configuration used to design a fuzzy logic controller using an evolutionary algorithm.}
\label{fig:gfsys}
\end{figure}
With all the chromosomes assigned a fitness value, genetic operators yield a new generation and the entire procedure is repeated until a stopping criterion, usually the number of iterations or the fitness value, is satisfied.



\subsection{Vehicle Model and Control Objective}
We use a ``bicycle model'' of a front-wheel-steering vehicle.\index{front-wheel-steering vehicle!bicycle model of}\index{bicycle model!of front-wheel-steering vehicle} Modeling of a ground vehicle is covered in detail in Subsection~\ref{subs:vehmod}. The standard coordinate system of the Society of Automotive Engineers~(SAE)\index{coordinate system!of the Society of Automotive Engineers}\index{Society of Automotive Engineers (SAE)}\index{SAE|see{Society of Automotive Engineers}} is used. A schematic of the model used is shown in Figure~\ref{fig:slc_model}.
\begin{figure}
\centering \epsfig{file=slc_model.eps,height=5.2cm}
\caption{Schematic of a bicycle model of a front-wheel-steering vehicle.}
\label{fig:slc_model}
\end{figure}
The modeling equations are:
\[
\left\{\begin{array}{lll}
{\dot x}_1 &=& -v_xx_3+\frac{1}{m}\left(F_{yf}+F_{yr}\right)\\
{\dot x}_2 &=& x_3\\
{\dot x}_3 &=&\frac{2}{I}(l_1F_{yf}-l_2F_{yr})\\
{\dot x}_4 &=& -v_xx_2-x_1,
\end{array}\right.
\]
where $x_1=\dot{y}$ is the lateral velocity of the vehicle,\index{lateral velocity!of vehicle} $x_2=\theta$ is the yaw angle, $x_3=\dot{\theta}$ is the yaw velocity,\index{yaw velocity} and $x_4=Y$ is the lateral position\index{lateral position} of the vehicle. Other parameters and their values are given in Table~\ref{tab:symbols}.
\begin{table}
\caption{Vehicle data and symbols used in modeling equations.}
\begin{center}
\begin{tabular}{cll}\hline
$\delta$ & steering angle of the front wheel & {}\\
$F_{yf}$ & lateral force of the front tire & {}\\
$F_{yr}$ &lateral force of the rear tire & {}\\
$m$ & total mass of the vehicle & 1280 kg \\
$I$ & moment of the inertia of the vehicle & 2500 kg$\cdot$m$^2$ \\
$l_1$ & distance from center of gravity to front axle & 1.2 m \\
$l_2$ & distance from center of gravity to rear axle & 1.22 m \\
$v_x$ & longitudinal velocity of the vehicle & 18.3 m/s \\\hline
\end{tabular}
\end{center}
\label{tab:symbols}
\end{table}
The lateral forces $F_{yf}=F_{yf}(\alpha_f)$ and $F_{yr}=F_{yr}(\alpha_r)$ are nonlinear functions of respective slip angles with a typical plot of a lateral force versus tire slip angle\index{tire slip angle} shown in Figure~\ref{fig:dugoffb}. Recall that the slip angles for the front and rear tires are
\[
\alpha_f=\delta-\frac{l_1x_3+x_1}{v_x}\quad\mb{and}\quad \alpha_r = \frac{l_2x_3-x_1}{v_x}.
\]
The input is the steering angle $\delta$ of the front wheel and the output is the lateral position of the vehicle, $x_4=Y$.

The step-lane-change maneuver\index{step-lane-change maneuver} along with cruising, traveling path selection, highway exiting and entering is one of the essential elements of automated vehicle control.\index{automated vehicle control} However, automation of these tasks requires a roadway reference system\index{roadway reference system} that involves sensing for lane detection.\index{lane detection} In our simulation, we assume that the road reference system is already installed. Our objective is to design a fuzzy logic controller that forces the vehicle to track a given reference lateral position signal $r$. The reference signal is generated by a reference signal generator\index{reference signal generator} whose transfer function is $1/(s^3+1.75s^2+2.15s+1)$. The tracking error $e$ is the difference between the reference lateral position signal $r$ and the actual lateral position $Y$ of the vehicle. The role of the fuzzy controller\index{fuzzy controller} is to produce the steering wheel angle $\delta$ using the information of the tracking error $e$ and the change in the tracking error $\dot e$ so that the vehicle tracks the reference signal. A SIMULINK block diagram of the closed-loop vehicle system\index{closed-loop vehicle system} is depicted in Figure~\ref{fig:sys}.
\begin{figure}
\centering \epsfig{file=sys.eps,width=15cm}
\caption{SIMULINK block diagram of the closed-loop vehicle system driven by a fuzzy logic controller.}
\label{fig:sys}
\end{figure}
In the following, we describe three different methods for an evolutionary algorithm~(EA) aided fuzzy logic controller~(FLC) design.\index{EA aided FLC design}

\subsection{Case 1: EA Tunes Fuzzy Rules}

In this case, we assume that the membership functions\index{membership functions} are fixed. In our simulation experiment, we use the membership functions shown in Figure~\ref{fig:mf_man}.
\begin{figure}
\centering \epsfig{file=mf_man.eps,height=8cm}
\caption{Fuzzy membership functions for $e$, $\dot{e}$, and $u$.}
\label{fig:mf_man}
\end{figure}
The labels of fuzzy sets\index{fuzzy sets} LN, N, Z, P, and LP denote large negative, negative, zero, positive, and large positive. There are five fuzzy sets for $e$ and five fuzzy sets for $\dot{e}$, hence we can have 25 different rules of the form:
\begin{center}
IF $e$ is LN AND $\dot e$ is LN then $\delta$ is SP\\ $\vdots$\\ IF $e$ is LP AND
$\dot e$ is LP then $\delta$ is SN.
\end{center}
These rules can be arranged into a rules table\index{rules table} as shown in
Figure~\ref{fig:plots_heu}(a). We mention here that the rules shown in
Figure~\ref{fig:plots_heu}(a) were arrived at heuristically. Using the rules table, we construct a FLC using the center average defuzzifier\index{center average defuzzifier} given by~(\ref{eq:product_de}).
\begin{figure}
\centering \epsfig{file=plots_heu.eps,height=12cm}
\caption{FLC constructed using a fuzzy rule base obtained heuristically: (a) fuzzy rules table (b) control surface (c) tracking performance.} \label{fig:plots_heu}
\end{figure}
The control surface\index{control surface} derived using the rules in Figure~\ref{fig:plots_heu}(a) is shown in Figure~\ref{fig:plots_heu}(b). Tracking performance for two
different initial conditions $Y(0)=5$~m and $Y(0)=-5$~m is illustrated
in Figure~\ref{fig:plots_heu}(c). We calculated tracking errors for each initial condition and labeled the tracking error for the first initial condition as $e_1$ while the tracking error of the second initial condition was labeled as $e_2$. We then constructed the performance index~($PI$)\index{performance index ($PI$)} as the sum of the integrals of the square of the error~(ISE)\index{square of the error (ISE)} over the simulation time,
\[
PI=\int^{50}_0 \left( e_1^2(t)+e_2^2(t)\right)dt.
\]
The value of the $PI$ in the above experiment is 107.04. The performance index, $PI$, was evaluated using the circuit shown in Figure~\ref{fig:evaluatepi}.
\begin{figure}
\setlength{\unitlength}{.03in}
\begin{center}
\begin{picture}(120,20)
\thicklines
\put(5,8){\em e} \put(10,10){\vector(1,0){15}}
\put(25,2){\framebox(20,16){$(\cdot)^2$}}
\put(45,10){\vector(1,0){15}}
\put(60,2){\framebox(15,16){\Large$\frac{1}{s}$}}
\put(75,10){\vector(1,0){15}} \put(93,8){$\int e^2 dt$}
\end{picture}
\end{center}
\caption{Circuit used to evaluate the performance index, $PI$.}
\label{fig:evaluatepi}
\end{figure}




Of course, we could include more error trajectories in the expression for the $PI$. However, for the sake of simplicity, we decided to use only two tracking error trajectories in our simulation experiments.

We next employ an EA to generate a fuzzy rules table automatically.\index{fuzzy rules!generated automatically} In order to represent the whole fuzzy rules table as a chromosome,\index{fuzzy rules table!as a chromosome} we encode the labels of the fuzzy singletons corresponding to the controller outputs as integers, where LN is represented by 1, N by 2, Z by 3, P by 4, and LP by 5. This encoding of the fuzzy singletons\index{fuzzy singletons} of the control signal $u$ is depicted in Figure~\ref{fig:mf_man}. The actual control signal, however, is produced as the weighted sum of the fuzzy singletons using the center average defuzzifier.\index{center average defuzzifier} Because there are 25 rules, each chromosome is a string composed of 25 integers. For example, the rules table of
Figure~\ref{fig:plots_heu}(a) can be represented as a chromosome of the form,
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline 5 & 5 & 5 & 4 & 3 & 5 & 5 & 4 & 3 & 2 & 4 & 4 & 3 & 2 & 2 &
4 & 3 & 2 & 1 & 1 & 3 & 2 & 1 & 1 & 1\\ \hline
\end{tabular}
\end{center}
In our experiment, we decided on the population size\index{population size} of 30. The roulette wheel method\index{roulette wheel method} was used to select chromosomes to the mating pool.\index{mating pool} We employed the single-point crossover\index{single-point crossover}\index{crossover!single-point} operator with the crossover rate\index{crossover rate} $p_c=0.9$. The role of the mutation operator\index{mutation operator} was to replace each element of a chromosome randomly with an integer from the set $\{1,\ldots, 5\}$. We used the mutation rate\index{mutation rate} $p_m=0.05$.

We also employed an elitist strategy\index{elitist strategy} where two best individuals of the current generation were passed unchanged to the next generation. The fitness function\index{fitness function} employed was the inverse of the $PI$ discussed above, that is,
\begin{equation}
Fitness = \frac{1}{PI}=\frac{1}{\int^{50}_0 \left(
e_1^2(t)+e_2^2(t)\right)dt}. \label{eqn:fitness_rule}
\end{equation}
Recall that $e_1(t)$ is the tracking error with the initial lateral position at 5 m, and $e_2(t)$ is the tracking error with the initial lateral position at $-5$~m. The input to the reference system was a square wave with the amplitude of 10~m and the frequency of 0.04~Hz. The rules table of the best FLC, after 50 iterations of the EA, is shown in Figure~\ref{fig:plots_rule}(a).
\begin{figure}
\centering \epsfig{file=plots_rule.eps,height=12cm}
\caption{The best FLC of the 50-th generation using a fuzzy rule base automatically generated by the EA: (a) fuzzy rules table (b) control surface (c) tracking performance.}
\label{fig:plots_rule}
\end{figure}
Figure~\ref{fig:plots_rule}(b) shows the fuzzy control surface\index{fuzzy control surface} based on the rules table of Figure~\ref{fig:plots_rule}(a). Two different initial conditions, $Y(0) = 5$~m and $Y(0) = -5$~m were
used to test the tracking performance of the fuzzy controller; Figure~\ref{fig:plots_rule}(c) depicts the results. The vehicle lateral position becomes almost indistinguishable from the reference position after about seven seconds. The value of the $PI$ is 97.10, which is an improvement over the performance of the FLC obtained using heuristically generated fuzzy rules.\index{fuzzy rules!generated heuristically}


\subsection{Case 2: EA Tunes Fuzzy Membership Functions}

In this case, fuzzy rules\index{fuzzy rules!fixed} are assumed to be known and fixed. The role of an EA is to tune the input and output membership functions.\index{membership functions!tuning using EA} We use the heuristically obtained rules table of Figure~\ref{fig:plots_heu}(a) as the given fuzzy rule base. The fuzzy sets for each input can be completely described using two real numbers if the corresponding fuzzy membership functions are symmetric with respect to the origin and, at each point of their universe of discourse, the sum of the values of the membership functions equals 1. The controller output fuzzy sets are also assumed to be symmetric with respect to the origin. Under the above assumptions, the controller output fuzzy singletons\index{fuzzy singletons} can be completely characterized using just two parameters. Hence, six real parameters can completely characterize all the input and output membership functions. So the chromosomes are strings of just six real parameters.\index{chromosomes!as strings of real parameters} For example, the fuzzy sets of Figure~\ref{fig:mf_man} can be represented as the string,
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}\hline
0.1 & 0.4 & 0.1 & 0.4 & 2 & 4 \\\hline
\end{tabular}
\end{center}
We proceed with a description the EA that operates on the above described fuzzy rule base.\index{rule base} The population size,\index{population size} the selection method,\index{selection method} the crossover operator\index{crossover operator} and its rate, the mutation operator\index{mutation operator} and its rate, and the fitness function, as well as other parameters, are the same as in the previous case. The resulting fuzzy sets, after 50 iterations of the EA, are shown
in Figure~\ref{fig:mf_mf}.
\begin{figure}
\centering \psfig{file=mf_mf.eps,height=8cm}
\caption{Input and output fuzzy membership functions of the best FLC of the 50-th generation obtained by tuning fuzzy membership functions.}
\label{fig:mf_mf}
\end{figure}
The control surface\index{control surface} and the tracking performance of the closed-loop system driven by the FLC constructed using fuzzy sets of Figure~\ref{fig:mf_mf} are depicted in Figure~\ref{fig:plots_mf}.
\begin{figure}
\centering \epsfig{file=plots_mf.eps,height=12cm}
\caption{The best FLC of the 50-th generation obtained by tuning fuzzy membership functions: (a) control surface (b) tracking performance.}
\label{fig:plots_mf}
\end{figure}
The $PI$ value is 43.51. The tracking error of the lateral position\index{lateral position!tracking error of} almost disappears in about four seconds.

\subsection{Case 3: EA Tunes Fuzzy Rules and Membership Functions}

We observed in the previous cases that changing fuzzy rules or modifying the shape of membership functions can improve the performance of the resulting FLC. So let us analyze an EA that simultaneously operates on fuzzy rules and fuzzy membership functions.\index{EA!simultaneously operating on fuzzy rules and membership functions} This approach is a modification of the technique described in~\cite{hwang99,kang00}.

In general, an $i$-th fuzzy rule has the form,\index{fuzzy rule}
\begin{center}
IF $x_1$ IS $F^i_{1}$ AND $\cdots$  AND $x_n$ IS $F^i_{n}$ THEN
$u_1$ IS $\theta^i_{1}$ AND $\cdots$  AND $u_m$ IS $\theta^i_{m}$
\end{center}
where $1\leq i\leq r$, $x_j$s are the state variables where $1\leq j \leq n$, and $u_k$s being the controller's outputs, where $1\leq k\leq m$. $F^i_j$s are fuzzy sets with trapezoidal membership functions, and $\theta^i_k$s are the controller's output variables whose membership functions are singletons. A membership
function, $\mu_{F^i_{j}}(x_j)$, of the fuzzy set, $F^i_j$, can be completely characterized using four parameters $c_{ij}$, $d_{ij}$, $l_{ij}$, and $r_{ij}$. It can be described as
\begin{equation}
\mu_{F^i_{j}}(x_j)=\left\{\begin{array}{lll}
\frac{x_j-l_{ij}}{c_{ij}-l_{ij}-d_{ij}}&\mb{if}& l_{ij} < x_j <
c_{ij}-d_{ij}\\ 1&\mb{if}& c_{ij}-d_{ij}\leq x_j <
c_{ij}+d_{ij}\\
-\frac{x_j-r_{ij}}{r_{ij}-c_{ij}-d_{ij}}&\mb{if}& c_{ij}+d_{ij}
< x_j < r_{ij}\\ 0&{}&\mb{otherwise,}
\end{array}\right.
\label{eq:trapez}
\end{equation}
or equivalently,
\begin{equation}
\mu_{F^i_{j}}(x_j)=\max\left(0,\min\left(1,\frac{x_j-l_{ij}}{c_{ij}-l_{ij}-d_{ij}},
-\frac{x_j-r_{ij}}{r_{ij}-c_{ij}-d_{ij}}\right)\right).
\label{eq:trapez2}
\end{equation}
Figure~\ref{fig:gfmod_mfs} shows a generic fuzzy membership
function of the form given by~(\ref{eq:trapez}) or, equivalently
by~(\ref{eq:trapez2}).
\begin{figure}
\centering \epsfig{file=gfmod_mfs.eps,height=3.5cm} \caption{An
example of a trapezoidal membership function.}
\label{fig:gfmod_mfs}
\end{figure}
The fuzzy singletons, on the other hand, are described as\index{fuzzy singleton!description of}
\begin{equation}
\mu_{\theta^i_{k}}(z)=\left\{\begin{array}{ll} 1&\mb{if } z =
\theta^i_{k}\\ 0&\mb{otherwise.}
\end{array}\right.
\label{eq:singleton}
\end{equation}
The controller's signal $u_k$ is calculated using the center
average defuzzification formula,
\begin{equation}
u_k=\frac{\sum^r_{i=1}w_i\theta^i_k}{\sum^r_{i=1}w_i},
\end{equation}
where $r$ is the number of fuzzy rules\index{number of fuzzy rules}\index{fuzzy rules!number of} and the firing strength,\index{firing strength} $w_i=\mu_{F^i_{1}}\times \mu_{F^i_{2}}\times\cdots
\times\mu_{F^i_{n}}$.

We present a method of representing a given fuzzy rule base as a chromosome.\index{fuzzy rule base!as a chromosome} In order to tune fuzzy rules and to alter the shape of membership functions, we require that a chromosome contains all the information pertaining to the structure of the fuzzy rules and the shape of membership functions explicitly. A chromosome is composed of two substructures: one substructure contains the fuzzy rules' structure matrix while the other substructure contains the information concerning the parameters of the membership functions of the fuzzy rule base. The fuzzy rules' structure matrix contains the information related to the absence or presence of the premise part of the fuzzy rule\index{premise part!of the fuzzy rule} corresponding to the particular controller input. Hence the elements of this matrix are from the set $\{0,1\}$.

The parameter matrix contains the parameters that describe fuzzy membership functions. As previously mentioned, there are four parameters for each fuzzy trapezoidal membership function and one for each fuzzy singleton. The parameter matrix has the same number of rows as the number of rules and the number of columns is equal four times the number of the controller inputs plus the number of the controller's outputs.
Figure~\ref{fig:rulepara}(a) shows an example of the rules' structure matrix and the parameter matrix for the controller with two inputs and one output whose operation is described by two IF-THEN rules. The fuzzy rules and membership functions corresponding to the rules' structure matrix and the parameter matrix in Figure~\ref{fig:rulepara}(a) are depicted in
Figure~\ref{fig:rulepara}(b). Note that a 0 in the second column of the first row of the rules' structure matrix means that the premise part of the first fuzzy rule corresponding to the second controller input is absent; that is, the first rule does not involve the controller's second input.
\begin{figure}
\centering \epsfig{file=rulepara_mat2.eps,height=8cm}
\caption{An example of a chromosome structure and its corresponding fuzzy rule base where the controller has two inputs, $x_1$ and $x_2$, and one
output $u$.}
\label{fig:rulepara}
\end{figure}


The fitness function\index{fitness function} has the form,
\begin{eqnarray}
Fitness &=& \frac{1}{PI+ \alpha\times(\mb{number of
rules})}\nonumber\\
&=&\frac{1}{\int^{50}_0 \left(e_1^2(t)+e_2^2(t)\right)dt+
\alpha\times(\mb{number of rules})},
\label{eqn:fitness}
\end{eqnarray}
where $e_1(t)$ is the tracking error with the initial lateral position at 5 m, and $e_2(t)$ is the tracking error with the initial lateral position at $-5$~m.  The reference signal used has the amplitude of 10~m and the frequency of 0.04~Hz over 50 seconds. The design parameter $\alpha$ is
employed to penalize for the excessive number of fuzzy rules. In our simulations, $\alpha=1$.

We next describe genetic operators\index{genetic operators} utilized by the EA. Three genetic operators are used with the probability $p_g$. We first discuss a genetic operator that modifies the number of fuzzy rules. Let $r_{new}$ be the number of fuzzy rules after the genetic operator is applied while $r_{old}$ be the number of fuzzy rules before the genetic operator is applied. Then,
\[
r_{new}=r_{old}+\Delta r,
\]
where $\Delta r$ is the correction term that we analyze next. The number of fuzzy rules cannot exceed a prespecified bound $r_{max}$. It is reasonable to assume that on the average $r_{old}$ is about $r_{max}/2$. Therefore, we wish to construct the correction term so that
\[
\Delta r \in [-r_{max}/2, r_{max}/2].
\]
Let $\Psi$ be a random variable\index{random variable!with standard normal distribution} with standard normal distribution $N(0,1)$. The probability that $\Psi$ will take its values in the interval $(-3,3)$ is 0.99. Let $\mb{integer}(\cdot)$ be the operator that produces the nearest integer value of its argument. Then, the random variable, $\mb{integer}(\Psi r_{max}/6)$, will take its value in the interval $( -\mb{integer}(r_{max}/2), \mb{integer}(r_{max}/2) )$ with the probability 0.99. In view of the above, we use the following expression for updating the number of fuzzy rules,
\begin{equation}
r_{new} = r_{old} +  \mb{integer}(\Psi r_{max}/6).
\label{eqn:rulemuta}
\end{equation}
This particular form of the genetic operator ensures that the variation range of the rules' number is reasonable with respect to the maximum number $r_{max}$ of the rules allowed. If the number of new rules, $r_{new}$, falls outside the allowed range, then we set $r_{new}=r_{old}$.

The second genetic operator acts upon the elements of the rules' structure matrix by taking the binary complement of each element, that is, $1\rightarrow 0$ and $0\rightarrow 1$ with probability $p_g$.

The third genetic operator acts upon the parameter matrix as per the following expression,
\begin{equation}
p_{ij}^{new} = p_{ij}^{old} + \eta_{ij}\Psi e^{1-F_{old}/F_{max}},
\label{eq:thirdgen}
\end{equation}
where $p_{ij}^{old}$ is the current value of the $(ij)$-th element of the parameter matrix, $p_{ij}^{new}$ is the updated value of the $(ij)$-th element of the parameter matrix, $\Psi$ is a random variable with $N(0,1)$, $F_{old}$ is the fitness value of the chromosome being operated upon, and $F_{max}$ is the maximum fitness value of the chromosomes of the current population. Each design parameter $\eta_{ij}$ was chosen according to the following recipe,
\[
\eta_{ij}=(\mb{maximal range of }p_{ij} - \mb{minimal range of }p_{ij})/20.
\]
The term $e^{1-F_{old}/F_{max}}$, in equation~(\ref{eq:thirdgen}),
allows a chromosome with low fitness to vary the elements of the parameter matrix in a larger range than for the chromosomes with high fitness values.

In our implementation of the EA, we employ the elitist strategy\index{elitist strategy} to ensure that  20\% of the fittest chromosomes in each generation are passed directly to the next generation without being acted upon by the genetic operators. The genetic operators, described above, are applied to the top 20\% of the chromosomes in the current generation to produce the remaining 80\% of the chromosomes of the new population.

After applying genetic operators to the chromosomes of the current generation, post-processing may be required to ensure that the resulting chromosomes are feasible. If a zero row in the rules' structure matrix is created, then that row is deleted and the corresponding row in the parameter matrix is also deleted. The following conditions must be satisfied by the parameter matrix for the chromosome containing this matrix to be feasible,
\[
\left.\begin{array}{lll}
d_{ij} &\ge& 0\\
l_{ij} &\leq& c_{ij}-d_{ij}\\
r_{ij} &\ge& c_{ij}+d_{ij}.
\end{array}\right\}
\]
If the above conditions are violated by a parameter, then that parameter is assigned a boundary value. For example, if $d_{ij} < 0$,
then $d_{ij}$ is set to 0, and if $l_{ij} > c_{ij}-d_{ij}$, then
$l_{ij}$ is set to $c_{ij}-d_{ij}$, etc.

The above described EA is used to construct a fuzzy logic controller\index{fuzzy logic controller!constructed using EA} for the front-wheel-steering ground vehicle model,\index{ground vehicle model!front-wheel-steering}\index{front-wheel-steering!ground vehicle model} whose dynamics are characterized by a certain type of symmetry that allows us to simplify computations. Specifically, suppose that in the rule base we have a fuzzy rule of the form,
\begin{center}
IF $x_1$ is $F_{1}$ and $x_2$ is $F_{2}$ THEN $u$ is $\theta$,
\end{center}
where the fuzzy sets $F_{1}$ and $F_{2}$ are specified by the membership functions $\mu_{F_1}(c_1,d_1,r_1,l_1)$ and $\mu_{F_2}(c_2,d_2,r_2,l_2)$,
respectively. The parameters $c_i$, $d_i$, $r_i$, and $l_i$, $i=1,2$, characterize membership functions of the form shown in Figure~\ref{fig:gfmod_mfs}. Because of the symmetry in the vehicle dynamics, the rule base also contains the rule
\begin{center}
IF $x_1$ is $\tilde{F}_{1}$ and $x_2$ is $\tilde{F}_{2}$ THEN $u$ is $-\theta$,
\end{center}
where $\tilde{F}_{1}$ and $\tilde{F}_{2}$ have membership functions
\[
\mu_{\tilde{F}_1}=\mu_{F_1}(-c_1,d_1,-l_1,-r_1)\quad\mb{and}\quad \mu_{\tilde{F}_2}=\mu_{F_2}(-c_2,d_2,-l_2,-r_2),
\]
respectively.


The size of the population in our simulations was 30, the maximum number of rules was $r_{max}=10$, and the probability $p_g$ was equal to $0.1$. A fuzzy rule base obtained after 100 iterations is shown in Figure~\ref{fig:gfmod_mf}. Only four rules were generated.
\begin{figure}[h]
\centering \epsfig{file=gfmod_mf.eps,height=8cm}
\caption{The best fuzzy rule base after 100 iterations of the EA operating on fuzzy membership functions and fuzzy rule base.}
\label{fig:gfmod_mf}
\end{figure}
The fuzzy rule base shown in Figure~\ref{fig:gfmod_mf} was then used to construct a fuzzy logic controller that was connected to the vehicle model in the closed-loop configuration shown in Figure~\ref{fig:sys}. The fuzzy control surface\index{control surface} constructed using the obtained fuzzy rule base is shown in Figure~\ref{fig:plots_gfmod}(a). The results of a typical simulation run, the fitness values of the best, average, and worst chromosomes in each generation are shown in
Figure~\ref{fig:plots_gfmod}(b), while the tracking performance is illustrated in Figure~\ref{fig:plots_gfmod}(c).
\begin{figure}[h]
\centering \epsfig{file=plots_gfmod.eps,height=12cm}
\caption{Summary of the simulation run of the closed-loop system driven by the fuzzy logic controller obtained from the fuzzy rule base shown in Figure~\ref{fig:gfmod_mf}: (a)~control surface (b)~fitness values of the best, average, and worst
chromosomes (c)~tracking performance.}
\label{fig:plots_gfmod}
\end{figure}
As in the previous cases, the simulations were performed for two different initial conditions, and excellent tracking was achieved in both cases. The value of the $PI$ was reduced to 38.06, which is lower than in the two previous designs.


\section{Notes}
Basic principles of genetics have their roots in the experiments performed by a nineteenth-century monk, Gregor Mendel\index{Mendel!basic principles of genetics} in his monastery garden in Brno now in the Czech Republic. Mendel worked with the garden pea ({\it Pisum sativum\/})\index{Pisum sativum!garden pea} because it possesses traits that are easily recognizable and appear in widely contrasting forms. In addition the garden pea lends itself well to controlled crosses. Mendel published his findings in 1866 in a paper titled, ``Experiments in Plant Hybridization''. According to Elseth and Baumgardner~\cite[page~3]{elseth84}, ``The publication of this paper marked the birth of genetics as a science. The principles stated in the paper still stand today as the cornerstone of modern genetics.'' The importance of Mendel's work went completely unrecognized for 34 years. ``In 1900, three botanists, Hugo de Vries (Dutch), Carl Correns (German), and Erich von Tschermak-Seysenegg (Austrian), simultaneously discovered Mendel's paper on heredity, recognized its importance, and cited it in their own publications.''~\cite[page~7]{elseth84}. We add that the term {\it gene\/}\index{gene} was coined by W.~L. Johannsen in 1903.

A seminal text on genetic algorithms by Holland~\cite{holland92} is definitely worth perusing. However, Holland's monograph is not intended as an introduction to genetic algorithms. A very well written introductory text to genetic algorithms is Mitchell's book~\cite{mitchell96}. A number of engineering design examples using genetic algorithms can be found in Gen and Cheng~\cite{gen97}. B\"{a}ck, Hammel, and Schwefel~\cite{back97} survey the history as well as the current state of evolutionary computation which is the term that encompasses evolution strategies, evolutionary programming, and genetic algorithms. Applications of genetic algorithms to automatic programming using evolving computer programs are described by Koza~\cite{koza92}. Automatic programming with genetic algorithms is called by Koza ``genetic programming''~(GP). Robust controller design with genetic algorithms are reported by Patton and Liu~\cite{patton94}, and Ge, Lee, and Zhu~\cite{ge96}. Kozek, Roska, and Chua~\cite{kozek93}, Pal and Bhandari~\cite{pal94}, and Zamparelli~\cite{zamparelli97} use genetic algorithms to train cellular neural networks. Kim and Myung~\cite{kim97} apply genetic algorithms to solve constrained optimization problems, while Koza et al.~\cite{koza97} employ genetic algorithms for the automatic design of electrical circuits. Yang, Hachino, and Tsuji~\cite{yang97} combine least-squares techniques with a genetic algorithm to on-line identification of time-delay systems. More recently, Loop et al.~\cite{loop10} use genetic algorithms to estimate regions of asymptotic stability of power electronics systems.

\section{Exercises}
\begin{exercise}\rm
Recall that the ``elitism'' is a genetic operator\index{elitism!a genetic operator} that forces the GA, or EA, to retain some number of the best individuals from the current population in the new population. Incorporate the elitism operator into your GA program and test it on the function~(\ref{eq:manminfu}). You may try, for example, to keep the best two individuals from the current population in the new population.
\end{exercise}
\begin{exercise}\rm
Implement the simple EA from Section~\ref{sec:simpleea} and test it on:
\begin{itemize}
\item[(i)] the function~(\ref{eq:manminfu});
\item[(ii)] the function $f(x)=x+|\sin(32 x)|$, where $0\le x<\pi$.
\end{itemize}
The function in~(ii) comes from Mitchell~\cite[page~9]{mitchell96}.
\end{exercise}
\begin{exercise}\rm
Implement the simple EA from Section~\ref{sec:simpleea} that maximizes the MATLAB's peaks function
\begin{eqnarray*}
f(x,y) &=&  3(1-x)^2\exp(-x^2 - (y+1)^2) - 10\left(\frac{x}{5} - x^3 - y^5\right)\exp(-x^2-y^2)\\
&{}&\quad - \frac{1}{3}\exp(-(x+1)^2 - y^2),
\end{eqnarray*}
where $-3\le x,y,\le 3$.
\end{exercise}
\begin{exercise}\rm
Implement the simple EA from Section~\ref{sec:simpleea} that maximizes the function
\[
f(x,y) =  3 + y\sin(0.7 x) + x\cos(1.5 y);
\]
where $-3.0\le x\le 12.1$ and $4.1\le y \le 5.8$.
\end{exercise}
\begin{exercise}\label{exer:tsp}\rm
Consider the classic traveling salesman problem~(TSP):\index{traveling salesman problem (TSP)}\index{TSP|see{traveling salesman problem}} Suppose a salesman must visit clients in different cities, and then return home. What is the shortest tour through those cities, visiting each one once and only once? Implement in MATLAB a simple variant of the genetic algorithm that optimizes the salesman route for 10 cities. How many different possible paths are there for 10 cities? Following Fogel~\cite[page~29]{fogel00}, use the following representation of any candidate solution: for a given map of the cities to be visited (including the salesman's home base), assign to each city a positive integer. Thus, for example, if there were six cities, one possible solution might be $\left[\begin{array}{cccccc}
1 & 2 & 3 & 4 & 5 & 6
\end{array}\right]$. The above represents an order of progression of the salesman's trip. Note that because the problem requires a round trip, the first coordinate of the vector representing the first city of the trip is also the last city visited by the salesman. Because the trip is a closed loop, it does not matter from which city the salesman starts his or her trip. As far as a crossover operator implementation is concerned, we recommend using Fogel's idea~\cite[page~29]{fogel00}, where a single-parent operator produces an offspring by inverting the visiting order between two randomly chosen coordinates of the parent vector. As an example of this implementation of the single-parent crossover operator, suppose that the parent is as above and that the randomly chosen coordinates, that is, the inversion points, are the second and fifth ones. The offspring would then be $\left[\begin{array}{cccccc}
1 & 5 & 3 & 4 & 2 & 6
\end{array}\right]$.

Plot a map with the cities and mark the obtained salesman's route. The coordinates of the cities are given in the following table.
\begin{table}[h]
\caption{Locations of the cities in Exercise~\ref{exer:tsp}.}
\begin{center}
\begin{tabular}{ccccccccccccccccccccc}\hline
$x$ &0.4306 &3.7094 &6.9330 &9.3582 &4.7758 &1.2910 &4.83831 & 9.4560 & 3.6774 &3.2849\\
$y$ & 7.7288 & 2.9727 &1.7785& 6.9080 &2.6394&4.5774 &8.43692 & 8.8150 &   7.0002& 7.5569
\\ \hline
\end{tabular}
\end{center}
%\label{tab:veh_p2}
\end{table}
\end{exercise}
\begin{exercise}\label{exer:flcga}\rm Consider a closed-loop system with a fuzzy logic controller~(FLC)\index{fuzzy logic controller (FLC)} in the loop shown in Figure~\ref{fig:flcga}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{flcga.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/flcga.bmp}
%\hskip1in\special{bmp:c:/mybook/flcga.bmp x=4in y=3in}
\caption{Closed-loop system with a fuzzy logic controller in the loop.}
\label{fig:flcga}
\end{figure}
Suppose that we use five membership functions for the FLC inputs and its output. A possible set of rules is given in Table~\ref{tab:rule_flcga}, where LN denotes large negative, N means negative, ZE---zero, P---positive, and LP refers to large positive.
\begin{table}
\caption{Rules matrix to be used to construct a fuzzy logic controller in Exercise~\ref{exer:flcga}.}
\begin{center}
\begin{tabular}{c||c|c|c|c|c} \hline
$\dot{e}$ & LN & N & ZE & P & LP \\ \cline{1-1}
$e$& {}&{}&{}&{}&{}\\ \hline \hline
LN & LN & LN & LN & N & ZE \\ \hline
N & LN & LN & N & ZE & P \\ \hline
ZE & N & N & ZE & P & P \\ \hline
P & N & ZE & P & LP & LP \\ \hline
LP & ZE & P & LP & LP & LP \\ \hline
\end{tabular}
\end{center}
\label{tab:rule_flcga}
\end{table}
Implement a variant of the simple EA from Section~\ref{sec:simpleea} that generates the set of rules to be used to synthesize the FLC. Assign to each fuzzy set a positive integer, for example, let LN be encoded as 1, N as 2, ZE as 3, P as 4, and LP as 5. Then, a candidate solution, that is shown in Table~\ref{tab:rule_flcga}, can be represented in a vector form as
\[
\left[\begin{array}{ccccccccccccccccccccccccc}
1& 1& 2& 2 & 3 &1&1&2&3&4&1&2&3&4&5&2&3&4&5&5&3&4&4&5&5
\end{array}\right],
\]
that is, a vector representation of the candidate solution given by the rules matrix is obtained using the stacking operation of the columns of the rules matrix.\index{candidate solution!stacked columns of the rules matrix} Use the one-point crossover\index{one-point crossover} operator similar to the one in the canonical GA. The mutation operation can be implemented by altering a randomly chosen coordinate in the given chromosome, that is, the vector representing a candidate solution. You can use the roulette-wheel method\index{roulette-wheel method} to determine the replacement. For example, if the seventh component was selected to be altered, then you can use the roulette-wheel to determine an integer from the set $\{1,2,3,4,5\}$ that will replace the ``1'' that occupies the seventh coordinate of the above chromosome. Chan, Lee, and Leung~\cite{kcc_chan97} used a similar version of the EA to generate fuzzy rules to synthesize a fuzzy logic target tracking algorithm.\index{fuzzy logic!target tracking algorithm}\index{target tracking algorithm!fuzzy logic based}
\end{exercise} 