**Scheduling in Manufacturing Systems – Ant Colony Approach**

Mieczysław Drabowski and Edward Wantuch

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/51487

## **1. Introduction**

[16] Csébfalvi, A. (2003). Optimal design of space structures with stability constraints. *Bontempi F (ed) System-Based Vision for Strategic and Creative Design, Vols 1-3: 2nd Inter‐ national Conference on Structural and Construction Engineering, September 23-26, Rome,*

[17] Csébfalvi, A. (2010). A Higher-Order Path-Following Method for Stability-Constrain‐ ed Optimization of Geometrically Nonlinear Shallow Trusses. O Allix, P Wriggers (ed) ECCM 2010, IV European Conference on Computational Mechanics, Palais des Congrès, Paris, France, May 16-21,2010: European Committee on Computational Sol‐

[18] Csébfalvi, A. An Improved ANGEL Algorithm for the Optimal Design of Shallow Truss Structures with Discrete Size and Continuous Shape Variables and Stability Constraints. *Erik Lund (ed) 9th World Congress on Structural and Multidisciplinary Opti‐*

[19] Hadi, M. N. S., & Alvani, K. S. (2003). Discrete Optimum Design of Geometrically Non-Linear Trusses using Genetic Algorithms, Seventh International Conference on The Application of Artificial Intelligence toCivil and Structural Engineering. *B.H.V.*

[20] Hanahara, K., & Tada, Y. Global Buckling Has to be Taken into Account for Optimal Design of Large-Scale Truss Structure. *Erik Lund (ed) 9th World Congress on Structural and Multidisciplinary Optimization: WCSMO-9. Shizuoka, Japan*, Shizuoka:paper ID:

[21] Lemonge, A. C. C., Barbosa, H. J. C., Fonseca, L. G., & Coutinho, A. L. G. A. (2010). A genetic algorithm for topology optimization of dome structures. *Helder C. Rodrigues (ed)2nd International Conference on Engineering Optimization, September 6-9, Lisbon, Por‐*

*mization: WCSMO-9. Shizuoka, Japan*, Shizuoka:paper ID: 100\_1, 1-8.

*Topping (Ed.), Civil-Comp Press, Stirling, Scotland*, paper 37.

*Italy, Leiden:Balkema Publishers*, 493-497, 9-05809-599-1.

ids, 2010.05.16-2010.05.21. , 1-7.

128 Ant Colony Optimization - Techniques and Applications

142\_1, 1-10.

*tugal*, paper ID: 01284, 1-15.

Scheduling problems, also in manufacturing systems [4], are described by following param‐ eters: the processing – computing – environments comprising processor (machines) set, oth‐ er resources comprising transportations and executions devices, processes (tasks) set and optimality criterion. We assume that processor set consists of *m* elements. Two classes of processors can be distinguished: dedicated (specialized) processors and parallel processors.

In production systems machines are regarded as dedicated rather than as parallel. In such a case we distinguish three types of dedicated processor systems: flow-shop, open-shop and job-shop. In the flow-shop all tasks have the same number of operations which are per‐ formed sequentially and require the same sets of processors. In the open-shop the order among the operations is immaterial. For the job-shop, the sequence of operations and the sets of required processors are defined for each process separately.

In the case of parallel processors each processor can execute any task. Hence, a task requires some number of arbitrary processors. As in deterministic scheduling theory [12] parallel processors are divided into three classes: identical processors – provided that all tasks are executed on all processors with that same productivity, uniform processors – if the produc‐ tivity depends on the processor and on the task, and unrelated processors – for which execu‐ tion speed depends on the processor and on the task. In each of the above cases productivity of the processor can be determined.

Apart from the processors the can be also a set of additional resources, each available in *m <sup>i</sup>* units.

The second parameter of the scheduling problem is the tasks system. The tasks correspond to the applications for manufactured goods. We assume that the set of tasks consists of n

© 2013 Drabowski and Wantuch; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Drabowski and Wantuch; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

tasks. For the whole tasks system it is possible to determine such feature as preemptability (or nonpreemptability) and existence (or nonexistence) of precedence constraints.

**2. Adaptation of ACO to solve the problems of scheduling**

The Ant Colony Optimization (ACO) algorithm [2] is a heuristics using the idea of agents (here: ants) imitating their real behavior. Basing on specific information (distance, amount of pheromone on the paths, etc.) ants evaluate the quality of paths and choose between them with some random probability (the better path quality, the higher probability it represents). Having walked the whole path from the source to destination, ants learn from each other by leaving a layer of pheromone on the path. Its amount depends on the quality of solution chosen by agent: the better solution, the bigger amount of pheromone is being left. The pheromone is then "vapouring" to enable the change of path chosen by ants and let them

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

131

ignore the worse (more distant from targets) paths, which they were walking earlier.

the division of tasks is important, but also their sequence.

**•** Vapouring factor of pheromone (from the range (0; 1));

ered and as a result, the best solution can long be uncovered.

**•** Number of agents (ants) in the colony;

mine relatively weak solutions.

doesn't have to be the best one).

The situation is similar for the vapouring factor:

**•** α – the amount of pheromone on the path;

been defined:

too big).

The result of such algorithm functioning is not only finding the solution. Very often it is the trace, which led us to this solution. It lets us analyze not only a single solution, but also per‐ mutations generating different solutions, but for our problems basing on the same division (i.e. tasks are scheduled in different order, although they are still allocated to the same pro‐ cessors). This kind of approach is used for solving the problems of synthesis, where not only

To adapt the ACO algorithm [24] to scheduling problems, the following parameters have

**•** For too big number of agents, the individual cycle of algorithm can last quite long, and the values saved in the table ("levels of pheromone") as a result of addition will deter‐

**•** On the other hand, when the number of agents is too small, most of paths will not be cov‐

**•** Too small value will cause that ants will quickly "forget" good solutions and as a result it can quickly come to so called *stagnation* (the algorithm will stop at one solution, which

**•** Too big value of this factor will make ants don't stop analyze "weak" solutions; further‐ more, the new solutions may not be pushed, if time, which has passed since the last solu‐ tion found will be long enough (it is the values of pheromone saved in the table will be

The ACO algorithm defines two more parameters, which let you balance between:

The process of choosing these parameters is important and should consider that:

Precedence constraints are represented as directed acyclic graphs (DAGs). Each task sepa‐ rately is described by a number of parameters. We enumerate tem in the following:


The optimality criteria constituting the third element of the scheduling problem are:


Due to the fact that scheduling problems and their optimizations are general NP-com‐ plete [10,25] we suggest meta-heuristic approach: Ant Colony Optimization and its com‐ parison with neural method and with polynomial algorithms for certain exemplary problems of task scheduling.

If a heuristic algorithm (such as ACO) finds an optimal solution to polynomial problems, it is probable that solutions found for NP-complete problems will also be optimal or least ap‐ proximated to optimal. ACO algorithm was tested with known polynomial algorithms and all of them achieved optimal solutions for those problems.

The comparisons utilized such polynomial algorithms as [3,5,12]:


For non-polynomial problems of tasks scheduling ACO algorithm was tested with list algo‐ rithms [12] (HLFET, HLFNET, SCFET, SCFNET), with PDF/HIS [18] for STG tasks and neu‐ ral approach [22].

## **2. Adaptation of ACO to solve the problems of scheduling**

The Ant Colony Optimization (ACO) algorithm [2] is a heuristics using the idea of agents (here: ants) imitating their real behavior. Basing on specific information (distance, amount of pheromone on the paths, etc.) ants evaluate the quality of paths and choose between them with some random probability (the better path quality, the higher probability it represents). Having walked the whole path from the source to destination, ants learn from each other by leaving a layer of pheromone on the path. Its amount depends on the quality of solution chosen by agent: the better solution, the bigger amount of pheromone is being left. The pheromone is then "vapouring" to enable the change of path chosen by ants and let them ignore the worse (more distant from targets) paths, which they were walking earlier.

The result of such algorithm functioning is not only finding the solution. Very often it is the trace, which led us to this solution. It lets us analyze not only a single solution, but also per‐ mutations generating different solutions, but for our problems basing on the same division (i.e. tasks are scheduled in different order, although they are still allocated to the same pro‐ cessors). This kind of approach is used for solving the problems of synthesis, where not only the division of tasks is important, but also their sequence.

To adapt the ACO algorithm [24] to scheduling problems, the following parameters have been defined:

**•** Number of agents (ants) in the colony;

tasks. For the whole tasks system it is possible to determine such feature as preemptability

Precedence constraints are represented as directed acyclic graphs (DAGs). Each task sepa‐

(or nonpreemptability) and existence (or nonexistence) of precedence constraints.

rately is described by a number of parameters. We enumerate tem in the following:

The optimality criteria constituting the third element of the scheduling problem are:

Due to the fact that scheduling problems and their optimizations are general NP-com‐ plete [10,25] we suggest meta-heuristic approach: Ant Colony Optimization and its com‐ parison with neural method and with polynomial algorithms for certain exemplary problems

If a heuristic algorithm (such as ACO) finds an optimal solution to polynomial problems, it is probable that solutions found for NP-complete problems will also be optimal or least ap‐ proximated to optimal. ACO algorithm was tested with known polynomial algorithms and

For non-polynomial problems of tasks scheduling ACO algorithm was tested with list algo‐ rithms [12] (HLFET, HLFNET, SCFET, SCFNET), with PDF/HIS [18] for STG tasks and neu‐

all of them achieved optimal solutions for those problems.

The comparisons utilized such polynomial algorithms as [3,5,12]:

**•** Number of operations,

130 Ant Colony Optimization - Techniques and Applications

**•** Resources requirements,

**•** Execution time,

**•** Ready time,

**•** Deadline,

**•** Weight.

**•** Schedule length,

**•** Mean flow time,

**•** Mean tardiness.

of task scheduling.

**•** Hu Algorithm,

**•** Baer Algorithm,

ral approach [22].

**•** Coffman-Graham Algorithm,

**•** Maximum lateness,

**•** Vapouring factor of pheromone (from the range (0; 1));

The process of choosing these parameters is important and should consider that:


The situation is similar for the vapouring factor:


The ACO algorithm defines two more parameters, which let you balance between:

**•** α – the amount of pheromone on the path;

**•** β – "quality" of the next step;

These parameters are chosen for specific task. This way, for parameters:


The process of agent involves:

ning sequence;

**4.** is it the last task?

**4.** release the task;

Idea of algorithm:

straints;

**7.** Go to 1;

*Example:*

**6.** If *G = Ø* END of algorithm;

*Algorithm:*

used:

the tasks to resources;

sources and running the tasks;

**1.** evaluation of current (incomplete) scheduling;

**3.** evaluation of the sequence obtained;

**5.** was it the last of available resources?

tasks, resources and times of tasks beginning.

**2.** allocation of task to the next of available resources;

**2.** With *S* select of tasks with the most strong of trace;

**5.** Update range of pheromone and remain of trace;

Two identical processors, digraph of seven tasks *Z <sup>i</sup> (t <sup>i</sup> )*, where *t <sup>i</sup>*

**1.** collecting information (from the tables of allocation) concerning allocation of tasks to re‐

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

133

**2.** drawing the next available task with the probability specified in the table of task run‐

**3.** drawing resources (processor) with the probability specified in the table of allocation

To evaluate the quality of allocation the task to processor, the following method is being

The calculative complexity of single agent is polynomial and depends on the number of

**1.** Construct *G* – structure of tasks non allocation and *S* – structure of tasks, which may be allocation in next step (for ex ample: begin: G = {Z1, Z2,…, Z7} and S = {Z1, Z2, Z3}); up‐

**3.** Allocate available of task as soon as possible and in accordance with precedence con‐

= time execution.

**4.** Remove selected of task with *G* and *S* and to add to list of tasks in memory of ant;

date range of pheromone and consideration of vapouring factor;


Having given the set of neighborhood *N* of the given point *i*, amount of pheromone on the path *τ* and the quality of passage from point *i* to point *j* as an element of the table *η* you can present the probability of passage from point *i* to *j* as [6,7]:

$$p\_{\boldsymbol{\eta}^{k}}^{(k)} = \begin{cases} \frac{\left[\boldsymbol{\pi}\_{\boldsymbol{\eta}}\right]^{\alpha} \left[\boldsymbol{\eta}\_{\boldsymbol{\eta}}\right]^{\beta}}{\sum\_{l \in N\_{l}^{k}} \left[\boldsymbol{\pi}\_{\boldsymbol{\eta}}\right]^{\alpha} \left[\boldsymbol{\eta}\_{\boldsymbol{\eta}}\right]^{\beta}} & \text{when } j \in N\_{l}^{k} \\\\ \frac{\left[\boldsymbol{\eta}\_{\boldsymbol{\eta}}\right]}{\left[\boldsymbol{\eta}\_{\boldsymbol{\eta}}\right]} & \text{else} \end{cases}$$

**Formula 1.** Evaluation of the quality of the next step in the ACO algorithm

In the approach presented here, the ACO algorithm uses agents to find three pieces of information:


Agents (ants) are searching for the solutions which are the collection resulting from the first two targets (they give the unique solution as a result). After scheduling, agents fill in two tables:


The process of agent involves:

**•** β – "quality" of the next step;

132 Ant Colony Optimization - Techniques and Applications

tion,

attendance),

information:

tables:

**•** the best sequence of tasks,

These parameters are chosen for specific task. This way, for parameters:

(ignorance of the level of pheromone on the path),

present the probability of passage from point *i* to *j* as [6,7]:

**Formula 1.** Evaluation of the quality of the next step in the ACO algorithm

**•** the best / the most beneficial division of tasks between processors,

**•** searching for the best possible solution for the given distribution.

**•** two-dimensional table representing allocation of task to the given processor,

**•** one-dimensional table representing the sequence of running the tasks.

of pheromone or the quality of solution.

**•** α > β there is bigger influence on the choice of path, which is more often exploited, **•** α < β there is bigger influence on the choice of path, which offers better solution,

**•** α = β there is balanced dependency between quality of the path and degree of its exploita‐

**•** α = 0 there is a heuristics based only on the quality of passage between consecutive points

**•** β = 0 there is a heuristics based only on the amount of pheromone (it is the factor of path

**•** α = β = 0 we'll get the algorithm making division evenly and independently of the amount

Having given the set of neighborhood *N* of the given point *i*, amount of pheromone on the path *τ* and the quality of passage from point *i* to point *j* as an element of the table *η* you can

In the approach presented here, the ACO algorithm uses agents to find three pieces of

Agents (ants) are searching for the solutions which are the collection resulting from the first two targets (they give the unique solution as a result). After scheduling, agents fill in two


To evaluate the quality of allocation the task to processor, the following method is being used:


The calculative complexity of single agent is polynomial and depends on the number of tasks, resources and times of tasks beginning.

Idea of algorithm:

*Algorithm:*


## *Example:*

Two identical processors, digraph of seven tasks *Z <sup>i</sup> (t <sup>i</sup> )*, where *t <sup>i</sup>* = time execution.

**•** Each component contains such number of neurons which equals the number of possible

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

135

**•** Assigning a specified value to a variable is the process of switching on a relevant neuron (neurons) and switching off the remaining ones in the component corresponding to this

**•** Constraints to the network are introduced by adding a negative weight connection be‐ tween neurons *('-1'*), symbolizing the variable values that cannot occur simultaneously.

Each neuron has its own table of connections and each connection contains its weight and the indicator for the connected neuron. A characteristic feature of the network is the diversi‐ ty of connections between neurons, but these never applied to all neurons [22,23]. It is a con‐ sequence of the fact that connections between neurons exist only when some constraints are imposed. The constraints existing in the discussed network model may be of the following

**•** Switching on a neuron means assigning the value *"1"* to its output.

**•** In the network there are additional neurons "the ones" that are switched on.

The method of constraints implementation shall be discussed upon examples [22].

Such net (Fig. 1.) blocks solution, in which Z1 = 1 as well as Z2 = 2 or Z3 = 3 as well as Z4 = 2.

Let us have two operations with unit execution times. The operation Z1 arrives at the system in time t = 1 and it is to be executed before the expiry of time t = 4. The operation Z2 arrives in time t = 1 and may be executed after the completion of operation Z1. A fragment of the net

*Neuron "one"* ('1') – a special neuron switched on permanently – is responsible for time con‐ straints. Introducing connections between such neuron and the relevant network neurons excludes a possibility of switching them on when searching for the solution. Task *Z <sup>1</sup>* cannot be scheduled in moment *0* and moment *4*, which corresponds to the assumption that this task arrives at the system at moment *1* and must be performed before moment *4*. Analogical process applies to operation *Z <sup>2</sup>*. The sequence constraints are executed by the connections between the network neurons. The figure shows (with dotted line) all the connections mak‐

**•** Switching off a neuron means assigning the *"0"* to its output.

for his case including all the connections is shown by Fig. 2.

ing the performance of task *Z <sup>2</sup>* before task *Z <sup>1</sup>* impossible.

values of each variable.

types: resources, time, order.

*Example 1:*

*Example 2:*

variable.

Parameters of ants' colony have been selected through experiments. Algorithm tuning is to select possibly best parameter values. This process demands many experiments which are conducted for different combinations of parameter values. For each combination of variable values, computation process has been repeated many times, and then an average result has been calculated. The same graphs type of STG, like at previous algorithms, have been ap‐ plied [18,20].

Selected algorithm parameters:


## **3. Adaptation of neural method to solve the problems of scheduling**

## **3.1. Neural network model**

The starting point for defining the neural network model for solving the problems of task scheduling and resource allocation are the assumptions for the constraint satisfaction prob‐ lem (CSP) [36,37]. CSP is the optimization problem which contains a certain set of varia‐ bles, sets of their possible values and constraints forced on the values of these variables [14,15]. On the basis of this problem assumption a network model of the following fea‐ tures is suggested:

**•** A neural network consists of components; each of them corresponds to another variable.


Each neuron has its own table of connections and each connection contains its weight and the indicator for the connected neuron. A characteristic feature of the network is the diversi‐ ty of connections between neurons, but these never applied to all neurons [22,23]. It is a con‐ sequence of the fact that connections between neurons exist only when some constraints are imposed. The constraints existing in the discussed network model may be of the following types: resources, time, order.

The method of constraints implementation shall be discussed upon examples [22].

## *Example 1:*

Parameters of ants' colony have been selected through experiments. Algorithm tuning is to select possibly best parameter values. This process demands many experiments which are conducted for different combinations of parameter values. For each combination of variable values, computation process has been repeated many times, and then an average result has been calculated. The same graphs type of STG, like at previous algorithms, have been ap‐

**•** a – number of ants; for number of tasks n < 50, a = 75 and for n>= 50, a = 1,5 x n

**3. Adaptation of neural method to solve the problems of scheduling**

The starting point for defining the neural network model for solving the problems of task scheduling and resource allocation are the assumptions for the constraint satisfaction prob‐ lem (CSP) [36,37]. CSP is the optimization problem which contains a certain set of varia‐ bles, sets of their possible values and constraints forced on the values of these variables [14,15]. On the basis of this problem assumption a network model of the following fea‐

**•** A neural network consists of components; each of them corresponds to another variable.

plied [18,20].

Selected algorithm parameters:

134 Ant Colony Optimization - Techniques and Applications

**3.1. Neural network model**

tures is suggested:

**•** γ – the pheromone evaporation coefficient = 0,08.

Such net (Fig. 1.) blocks solution, in which Z1 = 1 as well as Z2 = 2 or Z3 = 3 as well as Z4 = 2.

## *Example 2:*

Let us have two operations with unit execution times. The operation Z1 arrives at the system in time t = 1 and it is to be executed before the expiry of time t = 4. The operation Z2 arrives in time t = 1 and may be executed after the completion of operation Z1. A fragment of the net for his case including all the connections is shown by Fig. 2.

*Neuron "one"* ('1') – a special neuron switched on permanently – is responsible for time con‐ straints. Introducing connections between such neuron and the relevant network neurons excludes a possibility of switching them on when searching for the solution. Task *Z <sup>1</sup>* cannot be scheduled in moment *0* and moment *4*, which corresponds to the assumption that this task arrives at the system at moment *1* and must be performed before moment *4*. Analogical process applies to operation *Z <sup>2</sup>*. The sequence constraints are executed by the connections between the network neurons. The figure shows (with dotted line) all the connections mak‐ ing the performance of task *Z <sup>2</sup>* before task *Z <sup>1</sup>* impossible.

**3.** *Calculating the weighted sum of all neurons inputs.*

tions exist.

are switched *off*.

**4.** *Switching on the neuron with the highest input value.*

at the start, in each area three instances may happen:

switched on randomly, the remaining ones are switched *off*.

change and it remains switched *on*.

interrupted after a certain number of calls.

**5.** *Return to relaxation or – if there are no changes – exit from relaxation.*

**6.** *If there are connections (constraints) between the neurons that are switched on, each weight be‐*

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

137

The algorithm starts from allocating weight *'-1'* to all connections and then the start solution is generated. It is created by giving random values to the subsequent variables. This process takes place in a certain way: for each task i.e. in each area of the net such number of neurons is switched on as it is necessary for a certain task to be completed. The remaining, in the part which is responsible for its performance, neurons are being switched off. In the obtained re‐ sult there are many contradictions, specified by switching on the neurons where the connec‐

Therefore, the next step of the algorithm is the relaxation process, the objective of which is to "satisfy" the maximum numbers of limitations (backtracking). The objective is to obtain the result where the number of situations, where two switched on neurons of negative weight connection between them is the lowest. While switching on neurons with the biggest value

**•** If there is one neuron of the biggest value in the area, it is switched *on*; the remaining ones

**•** If there are more neurons, among which there is a previously switched one, there is no

**•** If there are more neurons, but there is no-one previously switched on, one of them is

A relaxation process finishes when the subsequent step does not bring any change and if all the requirements are met – the neurons between which a connection exist are not switched on – the right solution is found. If it is not still the case, it means that the algorithm found the local minimum and then the weight of each connection between two switched on neu‐ rons is decreased by *"1"* while its absolute value is being increased. It causes an increase in "interaction force' of this constraint which decreases the chance of switching *on* the same

After a certain number of iterations the network should consider all the constraints – provid‐ ing that there is the right solution, it should be found. Another factor is worth pointing out: in a relaxation process such an instance may occur where changes always happen. Then, this process might never be completed. Then a problem is solved in such a way that relaxation is

Search for a solution by algorithm consists of two stages. At the first one, which is described by the above presented algorithm, some activities are performed which lead to finding the right solution for the given specification. After finding such a solution, in consequence of

neurons in a relaxation process where we return in order to find the right solution.

*tween two switched on neurons is decreased by 1 and there is a return to relaxation.*

**Figure 1.** The example 1 of constraints.

**Figure 2.** The example 2 of constraints.

## **3.2. The algorithm description**

After entering the input data (the system specification), the algorithm constructs a neural network, the structure of which and the number of neurons composing it, depend upon the size and complexity of the instance of problem. We will name the part of the net allocated to this task – an area.

Constraints are introduced to the network by the execution of connections, occurring on‐ ly between the neurons corresponding to the values of variables which cannot occur simul‐ taneously.

The operation of the algorithm is the process of switching on appropriate neurons in each domain of network in order to satisfy the constraints imposed by the input data.

The algorithm course is as follows [36,38]:


**3.** *Calculating the weighted sum of all neurons inputs.*

**Figure 1.** The example 1 of constraints.

136 Ant Colony Optimization - Techniques and Applications

**Figure 2.** The example 2 of constraints.

**3.2. The algorithm description**

The algorithm course is as follows [36,38]:

**1.** *Allocating random values to consecutive variables.*

this task – an area.

**2.** *Network relaxation:*

taneously.

After entering the input data (the system specification), the algorithm constructs a neural network, the structure of which and the number of neurons composing it, depend upon the size and complexity of the instance of problem. We will name the part of the net allocated to

Constraints are introduced to the network by the execution of connections, occurring on‐ ly between the neurons corresponding to the values of variables which cannot occur simul‐

The operation of the algorithm is the process of switching on appropriate neurons in each

domain of network in order to satisfy the constraints imposed by the input data.


The algorithm starts from allocating weight *'-1'* to all connections and then the start solution is generated. It is created by giving random values to the subsequent variables. This process takes place in a certain way: for each task i.e. in each area of the net such number of neurons is switched on as it is necessary for a certain task to be completed. The remaining, in the part which is responsible for its performance, neurons are being switched off. In the obtained re‐ sult there are many contradictions, specified by switching on the neurons where the connec‐ tions exist.

Therefore, the next step of the algorithm is the relaxation process, the objective of which is to "satisfy" the maximum numbers of limitations (backtracking). The objective is to obtain the result where the number of situations, where two switched on neurons of negative weight connection between them is the lowest. While switching on neurons with the biggest value at the start, in each area three instances may happen:


A relaxation process finishes when the subsequent step does not bring any change and if all the requirements are met – the neurons between which a connection exist are not switched on – the right solution is found. If it is not still the case, it means that the algorithm found the local minimum and then the weight of each connection between two switched on neu‐ rons is decreased by *"1"* while its absolute value is being increased. It causes an increase in "interaction force' of this constraint which decreases the chance of switching *on* the same neurons in a relaxation process where we return in order to find the right solution.

After a certain number of iterations the network should consider all the constraints – provid‐ ing that there is the right solution, it should be found. Another factor is worth pointing out: in a relaxation process such an instance may occur where changes always happen. Then, this process might never be completed. Then a problem is solved in such a way that relaxation is interrupted after a certain number of calls.

Search for a solution by algorithm consists of two stages. At the first one, which is described by the above presented algorithm, some activities are performed which lead to finding the right solution for the given specification. After finding such a solution, in consequence of purpose function optimization there is a change of values for a certain criterion – in this case, decrease – then, the subsequent search for the right solution occur. In this case the search aims at a solution which possesses bigger constraints as the criteria value is sharper. Two criteria are taken into consideration for which a solution is being searched. It may be a cost function – where at the given time criterion, we search for the cheapest solution, or time function – where at the given cost criterion, we search for the quickest solution. Thus, the run of the algorithm is to seek a solution for smaller and smaller value of a selected criterion. However, if the algorithm cannot find the right solution for the recently modified criteria value of the algorithm, it returns to the previous criteria value for which it has found the right solution and modifies it by a smaller value.

For instance, if an algorithm has found the right solution for cost criterion which is e.g. 10, and it cannot find it for cost criteria which are 9, it tries to find a solution for cost 9.5 etc. In this way the program never finishes work, but all the time it tries to find a better solution in sense of a certain criterion. The user/designer of the system can interrupt its work at any moment if he/she considers the current solution given by an algorithm to be satisfying.

**Figure 3.** The task part for scheduling problems.

**Figure 4.** The resource part for scheduling problems.

Before selecting the quanta positions in the areas, algorithm has to calculate inputs for all the neurons. The neurons of the resource part are also connected to these inputs, as the number and the remaining places in recources have an impact on the setting which is going to "win" at a certain stage of computation. Thus, before an algorithm sets an exact task, it calculates the value of neuron inputs in resource part. The [r, i, k] neuron is switched on if at 'k' moment the resource 'r' is overloaded (too many tasks are using t), or it is not overload‐

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

139

The resource part:

In case of time criterion minimization, optimization goes at two planes. At the first one, sub‐ sequent neurons of the right side in task part of the network are connected to the neurons *"one"*, in this way fewer and fewer quanta is available for the algorithm of task scheduling which causes moving a critical line to the left and at the same time its diminishing. Howev‐ er, at the second, an individual quantum of time is being diminished; at each step an indi‐ vidual neuron will mean a smaller and smaller time passage.

## *The task part:*

Each area corresponds to one task (Fig. 3.). For further area, the best possible setting for the task is selected. Which setting 'wins' at the given stage and in the given area – this shall be determined by the sum of neuron outputs in the setting, i.e. the one that introduces the smaller number of contradictions. Moreover, it is checked if among the found set of the best solutions there is no previous one, then it is left.

A neuron at the *[i, k]* position corresponds to the presence of '*i*' task on the processor at the '*k*' moment. Between these neurons there are suitable inhibitory connections *(-1.0).*

If, for example, task 1 must be performed before task 2, for all the neuron pairs

*[1, k], [2, m]* there are inhibitory connections (denoting contradictions), if *k >= m* and if task 8 occurs in the system at moment 2, *"one"* neuron is permanently connected to neurons *[8, 0] and [8, 1]* (neuron which has 1.0 at the start which is permanently contradictory) and guaran‐ tees that in the final solution there is no quantum at moment 0 or 1.

We also take critical lines into account, which stand for time constraints that cannot be ex‐ ceeded by any allocated tasks – connecting 'one' will apply to the neurons of the right side of the network outside the critical line.

**Figure 3.** The task part for scheduling problems.

purpose function optimization there is a change of values for a certain criterion – in this case, decrease – then, the subsequent search for the right solution occur. In this case the search aims at a solution which possesses bigger constraints as the criteria value is sharper. Two criteria are taken into consideration for which a solution is being searched. It may be a cost function – where at the given time criterion, we search for the cheapest solution, or time function – where at the given cost criterion, we search for the quickest solution. Thus, the run of the algorithm is to seek a solution for smaller and smaller value of a selected criterion. However, if the algorithm cannot find the right solution for the recently modified criteria value of the algorithm, it returns to the previous criteria value for which it has found the

For instance, if an algorithm has found the right solution for cost criterion which is e.g. 10, and it cannot find it for cost criteria which are 9, it tries to find a solution for cost 9.5 etc. In this way the program never finishes work, but all the time it tries to find a better solution in sense of a certain criterion. The user/designer of the system can interrupt its work at any moment if he/she considers the current solution given by an algorithm to be satisfying.

In case of time criterion minimization, optimization goes at two planes. At the first one, sub‐ sequent neurons of the right side in task part of the network are connected to the neurons *"one"*, in this way fewer and fewer quanta is available for the algorithm of task scheduling which causes moving a critical line to the left and at the same time its diminishing. Howev‐ er, at the second, an individual quantum of time is being diminished; at each step an indi‐

Each area corresponds to one task (Fig. 3.). For further area, the best possible setting for the task is selected. Which setting 'wins' at the given stage and in the given area – this shall be determined by the sum of neuron outputs in the setting, i.e. the one that introduces the smaller number of contradictions. Moreover, it is checked if among the found set of the best

A neuron at the *[i, k]* position corresponds to the presence of '*i*' task on the processor at the

*[1, k], [2, m]* there are inhibitory connections (denoting contradictions), if *k >= m* and if task 8 occurs in the system at moment 2, *"one"* neuron is permanently connected to neurons *[8, 0] and [8, 1]* (neuron which has 1.0 at the start which is permanently contradictory) and guaran‐

We also take critical lines into account, which stand for time constraints that cannot be ex‐ ceeded by any allocated tasks – connecting 'one' will apply to the neurons of the right side

'*k*' moment. Between these neurons there are suitable inhibitory connections *(-1.0).*

If, for example, task 1 must be performed before task 2, for all the neuron pairs

tees that in the final solution there is no quantum at moment 0 or 1.

right solution and modifies it by a smaller value.

138 Ant Colony Optimization - Techniques and Applications

vidual neuron will mean a smaller and smaller time passage.

solutions there is no previous one, then it is left.

of the network outside the critical line.

*The task part:*

**Figure 4.** The resource part for scheduling problems.

The resource part:

Before selecting the quanta positions in the areas, algorithm has to calculate inputs for all the neurons. The neurons of the resource part are also connected to these inputs, as the number and the remaining places in recources have an impact on the setting which is going to "win" at a certain stage of computation. Thus, before an algorithm sets an exact task, it calculates the value of neuron inputs in resource part. The [r, i, k] neuron is switched on if at 'k' moment the resource 'r' is overloaded (too many tasks are using t), or it is not overload‐ ed, but setting the task of part 'i' at the moment defined by *'k'* would result in overloading. Neurons in resource part ( Fig. 4.) respond by their possible connection, resource overload‐ ed, if part of the task were set and at moment *'k';* therefore, neurons of resource part are con‐ nected to task inputs.

**4.** Allocation AAAA: there is no place at quantum 4 –algorithm adds the third processor

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

141

and allocates:

DDDDDD-----------





DDDDDD---EEEEE---

DDDDDD---EEEEE---

DDDDDD---EEEEE---



P1:DDDDDD---EEEEE---

may be estimated as follows:

P2:-C-BBB----AA-BB--

P3: ----AAAA---------

i – Number of tasks.

p – Number of processors.

Where:



**5.** Allocation EEEEE: there is place on the first processor :

**6.** Allocation AA there is place on the second processor

**7.** Allocation BB: there is place on the second processor:

The result of the operations on the processors is as follows:

Computational complexity of neural algorithm for task scheduling

An algorithm gives the right solution for the problems of known multi-nominal algorithms and also may be used for problems NP-complete. The complexity of one computation step

i \* 1( + ++ + + p \* k k\* 1( k \* p) \* m k \* r \* i \* p) p (2)

When in the resource part the neuron *" 'r' resource overload'* is switched *on*, as task *'i'* is set at moment *'k'* ", its signal *(1.0)* is transferred by weight (-1.0) to the neuron existing in the task part, which causes the negative input impulse (-1.0 \* 1.0) at the input which results in a con‐ tradiction.

In other words – it "disturbs" function *'compute\_in'* to set the task and at moment "k". Thus, in the network there are subsequent illegal situations implemented (constraints).

Each neuron *[r, i, k]* of the resource part is connected with neuron *[i, k]* from the task part, so a possibility of task existence at a given moment with concurrent resource overloading is 'in‐ hibited'.

Example:

Let us assume that there are five operations *A, B, C, D, E.*

Task part works as follows (a letter means a neuron switched on, sign '-' means a switched off neuron):


These operations should be allocated to a certain number of processors, so that one only op‐ eration would be performed on one processor at an exact moment:

**1.** Algorithm allocates ( at moment 0) fragment DDDDDD, adds a new processor ( the first) and allocates on it:

DDDDDD-----------

**2.** Allocation -C: for this moment (1) there is no place on the first processor, so algorithm adds the next processor and allocates an operation:

DDDDDD-----------


**3.** Allocation BBB: there is place on the second processor:

DDDDDD-----------


**4.** Allocation AAAA: there is no place at quantum 4 –algorithm adds the third processor and allocates:

DDDDDD-----------


ed, but setting the task of part 'i' at the moment defined by *'k'* would result in overloading. Neurons in resource part ( Fig. 4.) respond by their possible connection, resource overload‐ ed, if part of the task were set and at moment *'k';* therefore, neurons of resource part are con‐

When in the resource part the neuron *" 'r' resource overload'* is switched *on*, as task *'i'* is set at moment *'k'* ", its signal *(1.0)* is transferred by weight (-1.0) to the neuron existing in the task part, which causes the negative input impulse (-1.0 \* 1.0) at the input which results in a con‐

In other words – it "disturbs" function *'compute\_in'* to set the task and at moment "k". Thus,

Each neuron *[r, i, k]* of the resource part is connected with neuron *[i, k]* from the task part, so a possibility of task existence at a given moment with concurrent resource overloading is 'in‐

Task part works as follows (a letter means a neuron switched on, sign '-' means a switched

These operations should be allocated to a certain number of processors, so that one only op‐

**1.** Algorithm allocates ( at moment 0) fragment DDDDDD, adds a new processor ( the

**2.** Allocation -C: for this moment (1) there is no place on the first processor, so algorithm

in the network there are subsequent illegal situations implemented (constraints).

Let us assume that there are five operations *A, B, C, D, E.*

eration would be performed on one processor at an exact moment:

adds the next processor and allocates an operation:

**3.** Allocation BBB: there is place on the second processor:

nected to task inputs.

140 Ant Colony Optimization - Techniques and Applications

tradiction.

hibited'. Example:

off neuron):

first) and allocates on it:

DDDDDD-----------

DDDDDD-----------

DDDDDD-----------




**5.** Allocation EEEEE: there is place on the first processor :

DDDDDD---EEEEE---



**6.** Allocation AA there is place on the second processor

DDDDDD---EEEEE---



**7.** Allocation BB: there is place on the second processor:

DDDDDD---EEEEE---



The result of the operations on the processors is as follows:

P1:DDDDDD---EEEEE---

P2:-C-BBB----AA-BB--

P3: ----AAAA---------

Computational complexity of neural algorithm for task scheduling

An algorithm gives the right solution for the problems of known multi-nominal algorithms and also may be used for problems NP-complete. The complexity of one computation step may be estimated as follows:

$$\text{i } \begin{pmatrix} 1+\mathbf{p} \ \ast \ \mathbf{k}+\mathbf{k} \ \ast \end{pmatrix} \begin{pmatrix} 1+\mathbf{k} \ \ast \ \mathbf{p} \end{pmatrix} \text{ \begin{matrix} 1\\1 \end{matrix}} \begin{pmatrix} \ast \ \mathbf{p} \end{matrix} + \begin{matrix} \mathbf{p} \end{matrix} \tag{2}$$

Where:


k – Number of time quanta.

m – Number of all consecutive depend abilities between tasks.

r – Number of resources.

The largest complexity is generated by the process of increasing the number of tasks and an increase in the number of time quanta. Also, maximum number of processors and number of constitutive depend abilities in the introduced graph have a powerful effect on computa‐ tion. It is a pessimistic estimation; in practice, real complexity may be slightly smaller, but proportional to that. An algorithm itself is convergent i.e. step by step generates better and better solutions.

## **4. Tests of task scheduling algorithms**

## **4.1. The comparison with polynomial algorithms**

To show convergence of ACO algorithm towards optimum, one can compare their results with optimal results of already existing, precise, polynomial algorithms for certain exempla‐ ry problems of task scheduling. If a heuristic algorithm finds an optimal solution to polyno‐ mial problems, it is probable that solutions found for NP-complete problems will also be optimal or at least approximated to optimal. Heuristic algorithm described herein was test‐ ed with known polynomial algorithms and all of them achieved optimal solutions for those problems. The comparisons utilized such polynomial algorithms as:

**Figure 5.** Graph of tasks used for the comparison of ACO algorithm with Coffman and Graham algorithm (test prob‐

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

143

**•** Optimal scheduling for two processors obtained as a result of Coffman and Graham algo‐

**•** Optimal scheduling for two processors obtained as a result of ACO algorithm use (1a).

**•** Scheduling for three processors obtained as a result of Coffman and Graham algorithm

**Figure 6.** Optimal scheduling for two processors - Coffman and Graham algorithm (1a).

**Figure 7.** Optimal scheduling for two processors - ACO algorithm (1a).

lem no 1).

rithm use (1a).

use (1b).


Comparisons of ACO solutions with selected precise polynomial algorithms will be present‐ ed as an example.

Coffman and Graham algorithm

Scheduling of tasks which constitute a discretionary graph with singular performance times on two identical processors in order to minimize C*max.* Calculation complexity of the algo‐ rithm is O(n*<sup>2</sup>* ).

Test problem no 1:


k – Number of time quanta.

142 Ant Colony Optimization - Techniques and Applications

r – Number of resources.

better solutions.

m – Number of all consecutive depend abilities between tasks.

**4. Tests of task scheduling algorithms**

**4.1. The comparison with polynomial algorithms**

**•** Coffman – Graham Algorithm,

Coffman and Graham algorithm

**•** 2 identical processors (a), 3 identical processors (b).

**•** 15 tasks with singular performance times.

).

**•** Hu Algorithm,

**•** Baer Algorithm,

ed as an example.

rithm is O(n*<sup>2</sup>*

Test problem no 1:

**•** Graph with tasks:

problems. The comparisons utilized such polynomial algorithms as:

The largest complexity is generated by the process of increasing the number of tasks and an increase in the number of time quanta. Also, maximum number of processors and number of constitutive depend abilities in the introduced graph have a powerful effect on computa‐ tion. It is a pessimistic estimation; in practice, real complexity may be slightly smaller, but proportional to that. An algorithm itself is convergent i.e. step by step generates better and

To show convergence of ACO algorithm towards optimum, one can compare their results with optimal results of already existing, precise, polynomial algorithms for certain exempla‐ ry problems of task scheduling. If a heuristic algorithm finds an optimal solution to polyno‐ mial problems, it is probable that solutions found for NP-complete problems will also be optimal or at least approximated to optimal. Heuristic algorithm described herein was test‐ ed with known polynomial algorithms and all of them achieved optimal solutions for those

Comparisons of ACO solutions with selected precise polynomial algorithms will be present‐

Scheduling of tasks which constitute a discretionary graph with singular performance times on two identical processors in order to minimize C*max.* Calculation complexity of the algo‐

**Figure 5.** Graph of tasks used for the comparison of ACO algorithm with Coffman and Graham algorithm (test prob‐ lem no 1).

**•** Optimal scheduling for two processors obtained as a result of Coffman and Graham algo‐ rithm use (1a).


**Figure 6.** Optimal scheduling for two processors - Coffman and Graham algorithm (1a).

**•** Optimal scheduling for two processors obtained as a result of ACO algorithm use (1a).


**Figure 7.** Optimal scheduling for two processors - ACO algorithm (1a).

**•** Scheduling for three processors obtained as a result of Coffman and Graham algorithm use (1b).


**•** Optimal scheduling for two processors obtained as a result of Coffman and Graham algo‐

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

145

**•** Optimal scheduling for two processors obtained as a result of ACO algorithm use (2a).

**•** Non-optimal scheduling for three processors obtained as a result of Coffman and Graham

**•** Optimal scheduling for three processors obtained as a result of ACO algorithm use (2b).

For the problem of two processors (2a) both algorithms obtained optimal scheduling. In the case of three processors (2b) the Coffman and Graham algorithm did not find optimal sched‐

**Figure 11.** Optimal scheduling for 2 processors - Coffman and Graham algorithm (2a).

**Figure 13.** Non-optimal scheduling for 3 processors – Coffman and Graham algorithm (2b).

**Figure 12.** Optimal scheduling for 2 processors - ACO algorithm (2a).

**Figure 14.** Optimal scheduling for 3 processors - ACO algorithm (2b).

uling, whereas the ACO algorithm did find it without any difficulty.

rithm use (2a).

algorithm use (2b).

**Figure 8.** Problem scheduling for 3 processors - Coffman and Graham algorithm use (1b).

**•** Scheduling for three processors obtained as a result of ACO algorithm use (1b).


**Figure 9.** Optimal problem scheduling for 3 processors – ACO algorithm (1b).

For two processors (1a) ACO algorithm identical to Coffman and Graham algorithm ob‐ tained optimal scheduling. It was the same in the case of three processors (1b) – both algo‐ rithms obtained the same scheduling. Coffman and Graham algorithm is optimal only for two identical processors. For task graph under research it also found optimal scheduling for 3 identical processors.

Another test problem is shown by the non-optimality of Coffman and Graham algorithm for processor number greater than 2.

*Test problem no 2*:


**Figure 10.** Graph of tasks used for the comparison of ACO algorithm with Coffman and Graham algorithm (test prob‐ lem no 2).

**•** Optimal scheduling for two processors obtained as a result of Coffman and Graham algo‐ rithm use (2a).


**Figure 11.** Optimal scheduling for 2 processors - Coffman and Graham algorithm (2a).

**•** Optimal scheduling for two processors obtained as a result of ACO algorithm use (2a).


**Figure 12.** Optimal scheduling for 2 processors - ACO algorithm (2a).

**Figure 8.** Problem scheduling for 3 processors - Coffman and Graham algorithm use (1b).

**Figure 9.** Optimal problem scheduling for 3 processors – ACO algorithm (1b).

**•** 2 identical processors (a), 3 identical processors (b)

**•** 12 tasks with singular performance times.

3 identical processors.

*Test problem no 2*:

**•** Graph of tasks:

lem no 2).

processor number greater than 2.

144 Ant Colony Optimization - Techniques and Applications

**•** Scheduling for three processors obtained as a result of ACO algorithm use (1b).

For two processors (1a) ACO algorithm identical to Coffman and Graham algorithm ob‐ tained optimal scheduling. It was the same in the case of three processors (1b) – both algo‐ rithms obtained the same scheduling. Coffman and Graham algorithm is optimal only for two identical processors. For task graph under research it also found optimal scheduling for

Another test problem is shown by the non-optimality of Coffman and Graham algorithm for

**Figure 10.** Graph of tasks used for the comparison of ACO algorithm with Coffman and Graham algorithm (test prob‐

**•** Non-optimal scheduling for three processors obtained as a result of Coffman and Graham algorithm use (2b).


**Figure 13.** Non-optimal scheduling for 3 processors – Coffman and Graham algorithm (2b).

**•** Optimal scheduling for three processors obtained as a result of ACO algorithm use (2b).


**Figure 14.** Optimal scheduling for 3 processors - ACO algorithm (2b).

For the problem of two processors (2a) both algorithms obtained optimal scheduling. In the case of three processors (2b) the Coffman and Graham algorithm did not find optimal sched‐ uling, whereas the ACO algorithm did find it without any difficulty.

In another test example both algorithms were compared for the problem of task scheduling on two identical processors with singular and different performance times.

*Test problem no 3:*


**Figure 18.** Non-optimal problem scheduling for irregular task performance times – Coffman and Graham algorithm

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

147

**•** Optimal scheduling for irregular task performance times obtained as a result of ACO al‐

Both compared algorithms obtain optimal scheduling for the problem with regular (singu‐ lar) task performance times (3a). For different task performance times (3b) the Coffman and Graham algorithm does not obtain optimal scheduling, whereas the ACO algorithm

Scheduling of tasks with singular performance times which create a digraph of anti-tree

type on identical processors in order to minimize C*max*. Algorithm complexity is O(n).

**Figure 20.** Graph of tasks used for the comparison of ACO and Hu algorithms.

**Figure 19.** Optimal problem scheduling for irregular task performance times – ACO algorithm (3b)

(2b)

gorithm use (3b):

does obtain.

*Hu algorithm*

*Test problem no 1:*

**•** 3 identical processors,

**•** 11 tasks with singular performance times,

**•** Graph of tasks:

**Figure 15.** Graph of tasks used for the comparison of ACO algorithm with Coffman and Graham algorithm (test prob‐ lem no 3)

**•** Optimal scheduling for singular task performance times obtained as a result of Coffman and Graham algorithm use (3a).


**Figure 16.** Optimal problem scheduling for singular task performance times – Coffman and Graham algorithm (3a)

**•** Optimal scheduling for singular task performance times obtained as a result of ACO algo‐ rithm use (3a).


**Figure 17.** Optimal problem scheduling for singular task performance times – ACO algorithm (3a)

**•** Non-optimal scheduling for irregular task performance times obtained as a result of Coff‐ man and Graham algorithm use (3b).

**Figure 18.** Non-optimal problem scheduling for irregular task performance times – Coffman and Graham algorithm (2b)

**•** Optimal scheduling for irregular task performance times obtained as a result of ACO al‐ gorithm use (3b):


**Figure 19.** Optimal problem scheduling for irregular task performance times – ACO algorithm (3b)

Both compared algorithms obtain optimal scheduling for the problem with regular (singu‐ lar) task performance times (3a). For different task performance times (3b) the Coffman and Graham algorithm does not obtain optimal scheduling, whereas the ACO algorithm does obtain.

## *Hu algorithm*

In another test example both algorithms were compared for the problem of task scheduling

**•** 5 tasks with singular performance times (a), 5 tasks with different performance times (b)

**Figure 15.** Graph of tasks used for the comparison of ACO algorithm with Coffman and Graham algorithm (test prob‐

**•** Optimal scheduling for singular task performance times obtained as a result of Coffman

**Figure 16.** Optimal problem scheduling for singular task performance times – Coffman and Graham algorithm (3a)

**Figure 17.** Optimal problem scheduling for singular task performance times – ACO algorithm (3a)

**•** Optimal scheduling for singular task performance times obtained as a result of ACO algo‐

**•** Non-optimal scheduling for irregular task performance times obtained as a result of Coff‐

on two identical processors with singular and different performance times.

*Test problem no 3:*

**•** Graph of tasks:

lem no 3)

rithm use (3a).

and Graham algorithm use (3a).

man and Graham algorithm use (3b).

**•** 2 identical processors.

146 Ant Colony Optimization - Techniques and Applications

Scheduling of tasks with singular performance times which create a digraph of anti-tree type on identical processors in order to minimize C*max*. Algorithm complexity is O(n).

**Figure 20.** Graph of tasks used for the comparison of ACO and Hu algorithms.

*Test problem no 1:*



**Figure 24.** Optimal scheduling for problem 2 solved with Hu algorithm

**Figure 25.** Optimal scheduling for problem 2 solved with ACO algorithm

anti-tree type on two uniform processors in order to minimize Cmax.

**•** 2 uniform processors with speed coefficients b1 = 2, b2 =1.

**Figure 26.** Graph of tasks used for the comparison of ACO and Baer algorithms.

**•** 11 tasks with singular performance times.

uling obtained is optimal.

**•** Graph of tasks (anti-tree):

*Baer algorithm*

*Test problem:*

**•** Optimal scheduling for problem 2 obtained as a result of ACO algorithm use.

Both problems solved with Hu algorithm were also solved easily by ACO algorithm. Sched‐

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

149

Scheduling of indivisible tasks, with singular performance times, which create a graph of

**Figure 21.** Optimal scheduling for problem 1 solved with Hu algorithm.

**•** Optimal scheduling for problem 1 obtained as a result of ACO algorithm use.


**Figure 22.** Optimal scheduling for problem 1 solved with ACO algorithm.

*Test problem no 2:*


**Figure 23.** Graph of tasks used for the comparison of ACO and Hu algorithms (test problem no 2)

**•** Optimal scheduling for problem 2 obtained as a result of Hu algorithm use.


**Figure 24.** Optimal scheduling for problem 2 solved with Hu algorithm

**•** Optimal scheduling for problem 2 obtained as a result of ACO algorithm use.


**Figure 25.** Optimal scheduling for problem 2 solved with ACO algorithm

Both problems solved with Hu algorithm were also solved easily by ACO algorithm. Sched‐ uling obtained is optimal.

## *Baer algorithm*

**•** Graph of tasks (anti-tree):

148 Ant Colony Optimization - Techniques and Applications

*Test problem no 2:*

**•** 3 identical processors,

**•** Graph of tasks (anti-tree):

**•** Optimal scheduling for problem 1 obtained as a result of Hu algorithm use:

**•** Optimal scheduling for problem 1 obtained as a result of ACO algorithm use.

**Figure 23.** Graph of tasks used for the comparison of ACO and Hu algorithms (test problem no 2)

**•** Optimal scheduling for problem 2 obtained as a result of Hu algorithm use.

**Figure 21.** Optimal scheduling for problem 1 solved with Hu algorithm.

**Figure 22.** Optimal scheduling for problem 1 solved with ACO algorithm.

**•** 12 tasks with singular performance times,

Scheduling of indivisible tasks, with singular performance times, which create a graph of anti-tree type on two uniform processors in order to minimize Cmax.

*Test problem:*


**Figure 26.** Graph of tasks used for the comparison of ACO and Baer algorithms.

**•** Optimal scheduling for the problem solved with Baer algorithm, obtained as a result of ACO algorithm use.

In the chapter five types of list scheduling rules were compared: HLFET (Highest Levels First with Estimated Times), HLFNET (Highest Levels First with No Estimated Times), RANDOM, SCFET (Smallest Co-levels First with Estimated Times), SCFNET (Smallest Co-

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

151

The number of cases, in which the solution differs less than 5% from optimal solution, is ac‐ cepted as an evaluation criterion for the priority allocation rule. If for 90% of examined ex‐ amples the sub-optimal solution fit in the above range, the rule would be described as "almost optimal". This requirement is met only by HLFET rule, which gives results varying

**•** 12 tasks with different performance times: (Z0,1), (Z1,1), (Z2,7), (Z3,3), (Z4,1), (Z5,1),

levels First with No Estimated Times) [12].

(Z6,3), (Z7,2), (Z8,2), (Z9,1), (Z10,3), (Z11,1).

**Figure 28.** The graph of tasks used for the comparison of ACO and list algorithms

Scheduling obtained as a result of ACO algorithm operation.

**Figure 29.** Scheduling obtained with ACO algorithm.

from optimum by 4,4% on average.

**•** 2 identical processors.

**•** Graph of tasks:

Example:

**Figure 27.** Optimal scheduling for the problem solved with Baer algorithm, obtained as a result of ACO algorithm use

For the problem optimized with Baer algorithm, the ACO algorithm also obtains optimal solution.

## **4.2. Comparison of algorithms for non-polynomial problems of task scheduling**

## *4.2.1. NP- complete problem no 1:*

Scheduling *nonpreemptive*, independent tasks on identical processors for C*max* minimization.


**Table 1.** Scheduling nonpreemptive, independent tasks on identical processors.

For all problems under research algorithms found similar solutions. Only neural algorithm did worse – for the problem of scheduling 10 tasks on 3 identical processors, 20 tasks on 6 processors and 20 tasks on 8 processors as well ACO algorithm for the problem of schedul‐ ing 20 tasks on 3 identical processors.

## *4.2.2. NP-complete problem no 2:*

List scheduling with various methods of priority allocation

Because in general case the problem of scheduling dependent, nonpreemptable tasks is highly NP-complete, in some applications one can use polynomial approximate algorithms. Such algorithms are list algorithms.

In the chapter five types of list scheduling rules were compared: HLFET (Highest Levels First with Estimated Times), HLFNET (Highest Levels First with No Estimated Times), RANDOM, SCFET (Smallest Co-levels First with Estimated Times), SCFNET (Smallest Colevels First with No Estimated Times) [12].

The number of cases, in which the solution differs less than 5% from optimal solution, is ac‐ cepted as an evaluation criterion for the priority allocation rule. If for 90% of examined ex‐ amples the sub-optimal solution fit in the above range, the rule would be described as "almost optimal". This requirement is met only by HLFET rule, which gives results varying from optimum by 4,4% on average.

Example:

**•** Optimal scheduling for the problem solved with Baer algorithm, obtained as a result of

**Figure 27.** Optimal scheduling for the problem solved with Baer algorithm, obtained as a result of ACO algorithm use

For the problem optimized with Baer algorithm, the ACO algorithm also obtains optimal

Scheduling *nonpreemptive*, independent tasks on identical processors for C*max* minimization.

For all problems under research algorithms found similar solutions. Only neural algorithm did worse – for the problem of scheduling 10 tasks on 3 identical processors, 20 tasks on 6 processors and 20 tasks on 8 processors as well ACO algorithm for the problem of schedul‐

Because in general case the problem of scheduling dependent, nonpreemptable tasks is highly NP-complete, in some applications one can use polynomial approximate algorithms.

**Cmax Neural Algorithm** **Cmax ACO Algorithm**

**4.2. Comparison of algorithms for non-polynomial problems of task scheduling**

5 3 4 4 10 3 9 8 10 6 4 4 20 3 15 16 20 6 9 8 20 8 7 6

**Number of processors**

**Table 1.** Scheduling nonpreemptive, independent tasks on identical processors.

List scheduling with various methods of priority allocation

ACO algorithm use.

150 Ant Colony Optimization - Techniques and Applications

*4.2.1. NP- complete problem no 1:*

**Number of tasks**

ing 20 tasks on 3 identical processors.

Such algorithms are list algorithms.

*4.2.2. NP-complete problem no 2:*

solution.


**Figure 28.** The graph of tasks used for the comparison of ACO and list algorithms

Scheduling obtained as a result of ACO algorithm operation.

The length of obtained scheduling is compliant with the scheduling which was obtained by means of the best list scheduling available for this case and which is HLFET ("al‐ most optimal").

Algorithms were investigated by scheduling tasks represented with the same graph (50 STG

**Cmax Cmax** Number

**Cmax Cmax** Number

 2 267 267 388 267 412 4 155 157 4487 160 3339 8 155 154 89 155 112 16 155 155 10 155 8

In all researched problems algorithms under comparison found optimal solution. The only difference can be observed in the number of iterations needed to find an optimal solution.

For multiple criteria optimization in the following tests comparisons were made of compro‐ mise solutions for ACO algorithm with the results of neural algorithm. Optimization criteria were: time, cost and power consumption. Additional requirements and constraints were

adopted: maximum number of processors – 5, maximal cost – 3, maximal time – 25.

ACO algorithm needed less iterations than neural one to find the solution.

**5. Comparing ACO algorithm and neural algorithm**

 2 228 228 132 228 92 4 114 114 1401 114 925 8 57 61 4318 58 4442 16 48 48 58 48 33

**PDF/IHS Ant colony Neural**

of iterations

**PDF/IHS Ant colony Neural**

of iterations

**Cmax** Number

Scheduling in Manufacturing Systems – Ant Colony Approach

**Cmax** Number

of iterations

of iterations

http://dx.doi.org/10.5772/51487

153

tasks) on a different number of processors.

**Number of processors**

**Table 3.** Minimization of Cmax of dependent tasks (STG rand0008.stg)

**Table 4.** Minimization of Cmax of dependent tasks (STG rand0107.stg)

**Number of processors**

**Number of tasks**

**Number of tasks**

## *4.2.3. Comparison with PDF/HIS algorithm*

For research purposes a set of graphs was utilized from the website below: http:// www.kasahara.elec.waseda.ac.jp/schedule/index.html. Task graphs made available therein were divided into groups because of the number of tasks. Minimum scheduling length was calculated by means of PDF/HIS algorithm (Parallelized Depth First/ Implicit Heuristic Search) for every tasks graph. STG graphs are vectored, a-cyclic tasks graphs. Different task performance times, discretionary sequence constraints as well as random number of pro‐ cessors cause STG tasks scheduling problems to be NP-complete problems. Out of all solved problems heuristic algorithms under research did not find an optimal solution (assuming this is the solution obtained with PDF/IHS algorithm) only for three of them. However, re‐ sults obtained are satisfactory, because the deviation from optimum varies from 0,36% to 4,63% (table Tab 2).


**Table 2.** Comparison with PDF/IHS algorithm – the influence of tasks number

Algorithms were investigated by scheduling tasks represented with the same graph (50 STG tasks) on a different number of processors.


**Table 3.** Minimization of Cmax of dependent tasks (STG rand0008.stg)

The length of obtained scheduling is compliant with the scheduling which was obtained by means of the best list scheduling available for this case and which is HLFET ("al‐

For research purposes a set of graphs was utilized from the website below: http:// www.kasahara.elec.waseda.ac.jp/schedule/index.html. Task graphs made available therein were divided into groups because of the number of tasks. Minimum scheduling length was calculated by means of PDF/HIS algorithm (Parallelized Depth First/ Implicit Heuristic Search) for every tasks graph. STG graphs are vectored, a-cyclic tasks graphs. Different task performance times, discretionary sequence constraints as well as random number of pro‐ cessors cause STG tasks scheduling problems to be NP-complete problems. Out of all solved problems heuristic algorithms under research did not find an optimal solution (assuming this is the solution obtained with PDF/IHS algorithm) only for three of them. However, re‐ sults obtained are satisfactory, because the deviation from optimum varies from 0,36% to

most optimal").

4,63% (table Tab 2).

**STG Num-ber**

**of tasks** **Number of processors**

**Table 2.** Comparison with PDF/IHS algorithm – the influence of tasks number

**PDF/ IHS**

**Cmax Cmax** Number of

rand0008 50 2 281 281 117 0 281 80 0 rand0038 50 4 114 114 1401 0 114 818 0 rand0107 50 8 155 155 389 0 155 411 0 rand0174 50 16 131 131 180 0 131 190 0 rand0017 100 2 569 569 171 0 569 92 0 rand0066 100 4 253 253 4736 0 257 3644 1,58 rand0106 100 8 205 205 861 0 205 927 0 rand0174 100 16 162 162 265 0 162 216 0 rand0020 300 2 827 846 5130 2,30 830 4840 0,36 rand0095 300 8 382 394 5787 3,14 384 5253 0,52 rand0136 300 16 324 339 2620 4,63 324 3067 0

iterations

**Ant colony Neural**

**Cmax Number of itera-tons**

Difference [%]

Difference [%]

*4.2.3. Comparison with PDF/HIS algorithm*

152 Ant Colony Optimization - Techniques and Applications


**Table 4.** Minimization of Cmax of dependent tasks (STG rand0107.stg)

In all researched problems algorithms under comparison found optimal solution. The only difference can be observed in the number of iterations needed to find an optimal solution. ACO algorithm needed less iterations than neural one to find the solution.

## **5. Comparing ACO algorithm and neural algorithm**

For multiple criteria optimization in the following tests comparisons were made of compro‐ mise solutions for ACO algorithm with the results of neural algorithm. Optimization criteria were: time, cost and power consumption. Additional requirements and constraints were adopted: maximum number of processors – 5, maximal cost – 3, maximal time – 25.


**Chart 31.** Influence of number of tasks on time – minimization of time, cost and power consumption.

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

155

**Chart 32.** Influence of number of tasks on power consumption – minimization of time, cost and power consumption.

Additional requirements and constraints were adopted: maximum number of processors: 5,

**Chart 33.** Influence of number of tasks on cost – minimization of time, cost and power consumption with of cost

Results were illustrated in the following charts - Chart: 33, 34, and 35.

maximal cost: 8, maximal time: 50.

of memory.

**Table 5.** Comparison of Ant Colony and neural for minimization of time, cost and power consumption.

Results were illustrated on the following charts – Chart: 30, 31, and 32.

When comparing solutions obtained by the algorithms one cannot provide an unequivocal answer which of the optimization methods is better. Greater influence on the quality of of‐ fered solutions has the algorithm itself, especially its exploration capacity of admissible solu‐ tions space. When analyzing the graphs of interdependence between cost and task number, it appears that neural algorithm is more stable i.e. attempts to maintain low cost, despite an increase in the number of tasks. This results in worse task performance time what is very visible on the graph where time is contingent on the number of tasks. From power con‐ sumption analysis it is evident that ACO algorithm solutions are more beneficial.

**Chart 30.** Influence of number tasks on cost – minimization of time, cost and power consumption .

**Chart 31.** Influence of number of tasks on time – minimization of time, cost and power consumption.

**Number of tasks**

154 Ant Colony Optimization - Techniques and Applications

**Ant colony Neural Cost Time Power consumption Cost Time Power consumption**

5 1,75 6,75 9,26 1,00 3,90 4,39 10 1,50 6,20 35,47 1,50 8,50 11,61 15 2,75 18,00 22,96 2,00 16,00 17,85 20 1,75 12,83 35,45 2,00 22,50 20,31 25 2,00 14,50 51,25 2,00 22,00 28,93 30 2,75 16,90 63,58 2,50 23,00 35,01 35 2,00 18,00 78,30 2,50 24,67 36,12 40 2,75 17,75 104,68 2,50 17,00 72,52 45 2,25 21,75 99,50 2,50 18,67 79,02 50 2,25 23,88 113,26 2,50 21,00 88,57 55 2,50 25,00 164,58 2,50 22,50 95,33

**Table 5.** Comparison of Ant Colony and neural for minimization of time, cost and power consumption.

sumption analysis it is evident that ACO algorithm solutions are more beneficial.

**Chart 30.** Influence of number tasks on cost – minimization of time, cost and power consumption .

When comparing solutions obtained by the algorithms one cannot provide an unequivocal answer which of the optimization methods is better. Greater influence on the quality of of‐ fered solutions has the algorithm itself, especially its exploration capacity of admissible solu‐ tions space. When analyzing the graphs of interdependence between cost and task number, it appears that neural algorithm is more stable i.e. attempts to maintain low cost, despite an increase in the number of tasks. This results in worse task performance time what is very visible on the graph where time is contingent on the number of tasks. From power con‐

Results were illustrated on the following charts – Chart: 30, 31, and 32.

**Chart 32.** Influence of number of tasks on power consumption – minimization of time, cost and power consumption.

Additional requirements and constraints were adopted: maximum number of processors: 5, maximal cost: 8, maximal time: 50.

Results were illustrated in the following charts - Chart: 33, 34, and 35.

**Chart 33.** Influence of number of tasks on cost – minimization of time, cost and power consumption with of cost of memory.

**6. Conclusions**

**•**

for supporting planning process.

such anomalies. Take an example of this digraph of tasks:

**•** Diminishing of performance time for all the tasks ti

sults of different algorithms and of different runs.

than optimum scheduling (independently from choice list!):

Conducted research shows that presented algorithms for task scheduling obtain good solu‐ tions - irrespectively of investigated problem complexity. These solutions are considered op‐ timal or sub-optimal whose deviation from optimum does not exceed 5%. Heuristic algorithms proposed for task scheduling problems, especially ACO, should be a good tool

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

157

One should indicate a possible and significant impact of anomalies in task scheduling on the quality of the obtained results. The following examples [12] show a possibility of appearing

' = ti

For different problem instances, particular algorithms may achieve different successes; oth‐ ers may achieve worse results at different numbers of tasks. The best option is to obtain re‐

– 1 and the scheduling is longer

**Chart 34.** Influence of number of tasks on time – minimization of time, cost and power consumption with of cost of memory.


**Table 6.** Comparison of Ant colony and neural for minimization of time, cost and power consumption.

**Chart 35.** Influence of number of tasks on power consumption – minimization of time, cost and power consumption with memory cost

## **6. Conclusions**

**•**

**Chart 34.** Influence of number of tasks on time – minimization of time, cost and power consumption with of cost

10 6,50 2,00 37,99 4,50 6,00 7,52 20 1,50 18,50 33,15 4,00 11,00 19,07 30 5,90 23,00 82,41 5,00 14,00 30,98 40 7,00 23,00 121,56 5,00 18,00 37,33 50 4,25 16,20 186,05 5,00 21,00 49,99 60 2,50 32,00 175,24 5,00 25,00 60,20 70 2,50 38,00 167,59 5,00 29,00 69,35 80 3,25 37,00 183,67 5,00 32,00 79,19 90 4,25 28,60 328,73 5,00 36,00 98,39 100 6,75 30,33 336,36 5,50 39,00 101,62 110 4,25 41,80 435,77 5,00 43,00 115,53

**Table 6.** Comparison of Ant colony and neural for minimization of time, cost and power consumption.

**Chart 35.** Influence of number of tasks on power consumption – minimization of time, cost and power consumption

**Ant colony Neural Cost Time Power consumption Cost Time Power consumption**

of memory.

**Number of tasks**

156 Ant Colony Optimization - Techniques and Applications

with memory cost

Conducted research shows that presented algorithms for task scheduling obtain good solu‐ tions - irrespectively of investigated problem complexity. These solutions are considered op‐ timal or sub-optimal whose deviation from optimum does not exceed 5%. Heuristic algorithms proposed for task scheduling problems, especially ACO, should be a good tool for supporting planning process.

One should indicate a possible and significant impact of anomalies in task scheduling on the quality of the obtained results. The following examples [12] show a possibility of appearing such anomalies. Take an example of this digraph of tasks:


**•** Diminishing of performance time for all the tasks ti ' = ti – 1 and the scheduling is longer than optimum scheduling (independently from choice list!):

For different problem instances, particular algorithms may achieve different successes; oth‐ ers may achieve worse results at different numbers of tasks. The best option is to obtain re‐ sults of different algorithms and of different runs.

The goal of this scheduling is to find an optimum solution satisfying the requirements and constraints enforced by the given specification of the tasks and resources as well as criteria.

**Author details**

**References**

217-234.

*ley&Sons, Inc. New York*.

Mieczysław Drabowski1\* and Edward Wantuch1,2

1 Cracow University of Technology, Poland

\*Address all correspondence to: drabowski@pk.edu.pl

2 AGH University of Science and Technology, Poland

[1] Aggoune, R. (2004). Minimizing the makespan for the flow shop scheduling problem

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

159

[3] Błażewicz, J., Drabowski, M., & Węglarz, J. (1984). Scheduling independent 2-pro‐

[4] Błażewicz, J., Ecker, K., Pesch, E., Schmidt, G., & Węglarz, J. (1996). Scheduling Com‐

[5] Błażewicz, J., Ecker, K., Pesch, E., Schmidt, G., & Węglarz, J. (2007). Handbook on Scheduling, From Theory to Applications. *Springer-Verlag Berlin Heidelberg*.

[6] Blum, C. (2005). Beam-ACO- Hybridizing ant colony optimization with bean search:

[7] Blum, C., & Sampels, M. (2004). An ant colony optimization algorithm for shop scheduling problems. *Journal of Mathematical Modeling and Algorithm*, 3, 285-308.

[8] Breit, J., Schmidt, G., & Strusevich, V. A. (2003). Non-preemptive two-machine open shop scheduling with non-availability constraints. *Math. Method Opr. Res.*, 57(2),

[11] Cheng, T. C. E., & Liu, Z. (2003). Approximability of two-machine no-wait flowshop

[12] Coffman, E. G. Jr. (1976). Computer and Job-shop scheduling theory. *John Wi‐*

[13] Colak, S., & Agarwal, A. (2005). Non-greedy heuristiad augmented neural networks

An application to open shop schedling. *Comput. Oper. Res.*, 32, 1565-1591.

cessor tasks to minimize schedule length. *Inform. Proce. Lett.*, 18, 267-273.

with availability constraints. *Eur. J. Oper., Res.*, 153, 534-543.

puter and Manufacturing Processes. *Springer*.

[9] Brucker, P. (2004). Scheduling Algorithms. *Springer*.

[10] Brucker, P., & Knust, S. (2006). Complex Scheduling. *Springer*.

scheduling with availability constraints. *Opr. Res. Lett.*, 31, 319-322.

for the open-shop scheduling problem. *Naval Res. Logist.*, 52, 631-644.

[2] Ostfeld, Avi. (2011). Any Colony Optimization. *Rijeka, Croatia, InTech*.

As for the optimality criteria for the manufacturing system for better control, we shall as‐ sume its minimum cost, maximum operating speed and minimum power consumption.

We will apply multi-criteria optimization in sense of Pareto. The solution is optimized in sense of Pareto if it is not possible to find a better solution, regarding at least one criterion without deterioration in accordance to other criteria. The solution dominates other ones if all its features are better. Pareto ranking of the solution is the number of solutions in a pool which do not dominate it. The process of synthesis will produce a certain number of nondominated solutions. Although non-dominated solutions do not guarantee that they are an optimal Pareto set of solutions; nevertheless, in case of a set of suboptimal solutions, they constitute one form of higher order optimal set in sense of Pareto and they give, by the way, access to the problem shape of Pareto optimal set of solutions.

Let's assume that we want to optimize a solution of two contradictory requirements: the cost and power consumption Fig. 36.

While using a traditional way with one optimization function, it is necessary to contain two optimal criteria in one value. To do that, it is advisable to select properly the scales for the criteria; if the scales are selected wrongly, the obtained solution will not be optimal. The chart in the illustration shows where, using linearly weighed sum of costs, we will receive the solution which may be optimizes in terms of costs.

**Figure 36.** Set of optimal solutions in sense of Pareto.

Cost optimization, power and time consumption in the problem of scheduling is, undoubt‐ edly, the problem where the potential number of solutions in sense of Pareto is enormous.

Future research: others of instances of scheduling problems, and additional criteria, espe‐ cially in sense of Pareto and for dependable systems, are still open and this issue is now studied.

## **Author details**

The goal of this scheduling is to find an optimum solution satisfying the requirements and constraints enforced by the given specification of the tasks and resources as well as criteria. As for the optimality criteria for the manufacturing system for better control, we shall as‐ sume its minimum cost, maximum operating speed and minimum power consumption.

We will apply multi-criteria optimization in sense of Pareto. The solution is optimized in sense of Pareto if it is not possible to find a better solution, regarding at least one criterion without deterioration in accordance to other criteria. The solution dominates other ones if all its features are better. Pareto ranking of the solution is the number of solutions in a pool which do not dominate it. The process of synthesis will produce a certain number of nondominated solutions. Although non-dominated solutions do not guarantee that they are an optimal Pareto set of solutions; nevertheless, in case of a set of suboptimal solutions, they constitute one form of higher order optimal set in sense of Pareto and they give, by the way,

Let's assume that we want to optimize a solution of two contradictory requirements: the cost

While using a traditional way with one optimization function, it is necessary to contain two optimal criteria in one value. To do that, it is advisable to select properly the scales for the criteria; if the scales are selected wrongly, the obtained solution will not be optimal. The chart in the illustration shows where, using linearly weighed sum of costs, we will receive

Cost optimization, power and time consumption in the problem of scheduling is, undoubt‐ edly, the problem where the potential number of solutions in sense of Pareto is enormous.

Future research: others of instances of scheduling problems, and additional criteria, espe‐ cially in sense of Pareto and for dependable systems, are still open and this issue is now

access to the problem shape of Pareto optimal set of solutions.

the solution which may be optimizes in terms of costs.

**Figure 36.** Set of optimal solutions in sense of Pareto.

studied.

and power consumption Fig. 36.

158 Ant Colony Optimization - Techniques and Applications

Mieczysław Drabowski1\* and Edward Wantuch1,2


## **References**


[14] Dechter, R., & Pearl, J. (1988). Network-based heuristic for constraint satisfaction problems. *Artificial Intelligence*, 34, 1-38.

[29] Lee, C. Y. (2004). Machine scheduling with availably constraints. *Leung J.Y.T. Hand‐*

Scheduling in Manufacturing Systems – Ant Colony Approach

http://dx.doi.org/10.5772/51487

161

[30] Meseguer, P. (1989). Constraint satisfaction problems: An overview. *AICOM*, 2, 3-17.

[31] Montgomery, J., Fayad, C., & Petrovic, S. (2006). Solution representation for job shop

[32] Morton, T. E., & Pentico, D. W. (1993). Heuristic Scheduling System. *Wiley, New York*.

[33] Nuijten, W. P. M., & Aarts, E. H. L. (1996). A computational study of constraint satis‐

[34] Pinedo, M. (2001). Scheduling Theory, Algorithms, and Systems,. *Prentice Hall, Engle‐*

[35] Taillard, E. (1993). Benchmarks for basic scheduling problems. *European J. Oper. Res.*,

[37] Wang, C. J., & Tsang, E. P. K. (1991). Solving constraint satisfaction problems using neural-networks,. *IEEE Second International Conference on Artificial Neural Networks*.

[38] Xu, J., & Parnas, D. L. (1993). On Satisfying Timing Constraints in Hard-Real-Time

[36] Tsang, E. (1993). Foundations of Constraint Satisfaction. *Academic Press, Essex*.

Systems. *IEEE Trans. on Software Engineering*, 19(1), 70-84.

faction for multiple capacitated job shop scheduling. *European J. Oper. Res.*, 90,

scheduling problems in ant colony optimization. *LNCS*, 4150, 484-491.

*book of Scheduling, CRC Press*, 22, 1-22.

269-284.

*wood Cliffs, N.J.*

64, 278-285.


[29] Lee, C. Y. (2004). Machine scheduling with availably constraints. *Leung J.Y.T. Hand‐ book of Scheduling, CRC Press*, 22, 1-22.

[14] Dechter, R., & Pearl, J. (1988). Network-based heuristic for constraint satisfaction

[15] Dorndorf, U., Pesch, E., & Phan-Huy, T. (2000). Constraint propagation techniques

[16] Drabowski, M., & Wantuch, E. (2006). Coherent Concurrent Task Scheduling and Re‐ source Assignment in Dependable Computer Systems Design,. *International Journal of Reliability Quality and Safety Engineering, World Scientific Publishing,*, 13(1), 15-24.

[17] Drabowski, M. (2007). Coherent synthesis of heterogeneous system- an ant colony optimization approach. *Proceedings of Artificial Intelligence Studies, Vol.4 (89)/2007, sup‐*

[18] Drabowski, M. (2007). The ant colony in par-synthesis of computer system. *Proceed‐ ings of the 11th IASTED International Conference on Artificial Intelligence and Soft Com‐*

[19] Drabowski, M. (2007). Coherent synthesis of heterogeneous system- an ant colony

[20] Drabowski, M. (2007). An Ant Colony Optimization to scheduling tasks on a grid.

[21] Drabowski, M. (2008). Solving Resources Assignment and Tasks Scheduling Prob‐

[22] Drabowski, M. (2008). Neural networks in optimization scheduling resources and processes for management on a grid. *Polish Journal of Environmental Studies*, 17(4C).

[23] Drabowski, M. (2009). Ant Colony and Neural method for scheduling of complex of operations and resources frameworks- comparative remarks. *Proceedings of the IAST‐ ED International Conference on Computational Intelligence, Honolulu, USA, ACTA Press,*

[24] Drabowski, M. (2011). Ant Colony Optimization for coherent synthesis of computer system. *Ostfeld A., (ed.) Ant Colony Optimization, InTech, Croatia, Austria, India*,

[25] Garey, M. R., & Johnson, D. S. (1979). Computers and intractability: A guide to the

[26] Ha, S., & Lee, E. A. (1997). Compile-Time Scheduling of Dynamic Constructs in Data‐

[27] Leung, J. Y. T. (2004). Handbook on Scheduling: Algorithms, Models and Perform‐

[28] Lee, C. Y. (1996). Machine scheduling with an availably constraint. *J. Global Optim.*, 9,

*puting, Palma de Mallorca, ACTA Press, Anaheim, USA*, 244-249.

lems using Neural Networks. *Artificial Intelligence Studies*, 2.

theory of NP-completeness,. *San Francisco, Freeman*.

ance Analysis,. *Chapman&Hall, Boca Raton*.

flow Program Graphs,. *IEEE Trans. On Computers*, 46(7).

optimization approach. *Studia Informatica*, 2.

*Polish Journal of Environmental Studies*, 16(5B).

for disjunctive scheduling problems. *Artificial Intelligence*, 122, 189-240.

problems. *Artificial Intelligence*, 34, 1-38.

160 Ant Colony Optimization - Techniques and Applications

*ported by IEEE, Siedlce*, 65-74.

*Anaheim, USA*, 91-97.

179-204.

363-384.


**Chapter 7**

**Provisional chapter**

**Traffic-Congestion Forecasting Algorithm Based on**

**Traffic-Congestion Forecasting Algorithm Based on**

The growth of **i**ntelligent **t**ransport **s**ystems (ITS) has recently been quite fast and impressive, and various kinds of studies on ITS from the viewpoint of artificial intelligence have also been done [1][2][3][4][5]. However, there are still many problems that need to be solved and alleviating traffic congestion is one of the main issues. Reducing traffic congestion is quite urgent because the amount of money lost due to congestion within only 1 km in Tokyo has reached as much as 400 million yen per year. To alleviate this situation, two traffic-control systems called the "Vehicle **I**nformation and **C**ommunication System (VICS)" and "the probe

VICS is a telecommunication system that transmits information such as that on traffic congestion and the regulation of traffic by detecting car movements with sensors installed on the road [6]. Information on car movements and that on forecasting traffic congestion are analyzed at the VICS center in real time and then the information from the center is displayed on equipment, such as car-navigation systems installed in individual cars (see Fig. 1).

> ©2012 Kurihara, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Kurihara; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

© 2013 Kurihara, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

distribution, and reproduction in any medium, provided the original work is properly cited.

**Pheromone Communication Model**

**Pheromone Communication Model**

Additional information is available at the end of the chapter

**c**ar **s**ystem (PCS)" are currently in operation in Japan.

Additional information is available at the end of the chapter

Satoshi Kurihara

**1. Introduction**

**Figure 1.** VICS

Satoshi Kurihara

http://dx.doi.org/10.5772/52563

**Provisional chapter**

## **Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model Pheromone Communication Model**

**Traffic-Congestion Forecasting Algorithm Based on**

Satoshi Kurihara Additional information is available at the end of the chapter

Satoshi Kurihara

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/52563

## **1. Introduction**

The growth of **i**ntelligent **t**ransport **s**ystems (ITS) has recently been quite fast and impressive, and various kinds of studies on ITS from the viewpoint of artificial intelligence have also been done [1][2][3][4][5]. However, there are still many problems that need to be solved and alleviating traffic congestion is one of the main issues. Reducing traffic congestion is quite urgent because the amount of money lost due to congestion within only 1 km in Tokyo has reached as much as 400 million yen per year. To alleviate this situation, two traffic-control systems called the "Vehicle **I**nformation and **C**ommunication System (VICS)" and "the probe **c**ar **s**ystem (PCS)" are currently in operation in Japan.

VICS is a telecommunication system that transmits information such as that on traffic congestion and the regulation of traffic by detecting car movements with sensors installed on the road [6]. Information on car movements and that on forecasting traffic congestion are analyzed at the VICS center in real time and then the information from the center is displayed on equipment, such as car-navigation systems installed in individual cars (see Fig. 1).

**Figure 1.** VICS

Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Kurihara; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Kurihara, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

©2012 Kurihara, licensee InTech. This is an open access chapter distributed under the terms of the Creative

PCS also provides information on car movements and that on forecasting traffic congestion to individual drivers the same as VICS does. Different from VICS, this system collects traffic information from all cars, which are considered to be movable sensor units. Each car has a telecommunication unit and transmits several kinds of information such as position, velocity, and the status of the car to the central server. Then, the calculated car-movement and traffic-congestion-forecasting information are analyzed and the information from the center is displayed on equipment, such as car-navigation systems.

needs a certain amount of time. Therefore, we propose a new congestion-forecasting system that can react to dynamically changing traffic conditions based on a coordination mechanism using the pheromone-communication model. Its main feature is to be able to forecast short-term congestion one or two minutes ahead. There have been many studies on ITS

Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model

http://dx.doi.org/10.5772/52563

165

Section 2 discusses traffic-congestion information for car agents and proposes a method of forecasting congestion that uses a multi-agent coordination mechanism for road agents set up at intersections based on a pheromone-communications model, which adaptively responds to increasing amounts of congestion. Section 3 discusses our tests to verify the

Congestion-forecasting technology is one of the main elements of ITS. Up to now, several

• Long-term forecasting of congestion: A method of statistically analyzing past traffic data,

• Short-term forecasting of congestion: A method of forecasting congestion a few minutes

Although it can effectively make forecasts under regular-congestion conditions that have originated from car and road situations, a large amount of past data is necessary for analysis. Moreover, it has weaknesses in forecasting under irregular-congestion conditions, such as those experienced during the Golden-week holidays in Japan and the Christmas-holiday

VICS and PCS essentially belong to the second classification, and is excellent at short-term forecasting of congestion. Yet, in the current VICS and PCS that is operating in Japan, data from each car is collected at the central server and all calculations are done there. Consequently, it is difficult to supply real-time information due to bottlenecks and time

For solving these problems, a short-term system of forecasting congestion based on distributed processing is adequate. This paper focuses on "roads", and we propose a MAS consisting of many road agents. These road agents are set up at every intersection, and they coordinate locally with one another to forecast congestion. In this research, we adopted the

Pheromone communications are based on the behavior of social insects like ants and bees and are applied as a model that is used to adaptively respond to dynamic changes in the environment in various applications [12] (see Fig. 2). Our method is an effective way of forecasting congestion in which each road agent generates pheromone information on its own road unit and exchanges this information with its neighboring agents. In a related study, Ando et al. investigated the forecasting of congestion in a local area a short time after pheromones had evaporated and diffused [11]. However, all drivers need to have the same probe-car system installed in their vehicles. Even though there has been some discussion on

[8][9][10], but there have been few on the forecasting of short-term congestion.

basic effectiveness of this method. Finally, we conclude this paper in Section 4.

**2. Method of forecasting congestion based on pheromone-**

methods have been proposed, two of which are classified below.

and discovering a pattern where congestion has occurred [13].

lags in communicating information and the centralized calculations.

pheromone-communications model as the mechanism for coordination.

**communications model**

season in the U.S.

ahead by using real-time information.

Though these two systems are currently operated in Japan, the system structure of both systems is top down and centralized, so the reaction to dynamic changes in traffic congestion and occurrence of accidents is usually delayed and serious problems can occur when the central server is down. In other words, there is a lack of real-time features and of robustness in these systems.

On the other hand, traffic-control systems like ITS and PCS essentially have interesting features for the coordination mechanisms of multi-agent systems (MASs). The coordination mechanisms of MAS can generally be divided into two types: direct and indirect. In the former, precise coordination can be achieved, but when the number of agents becomes excessive the load of coordination becomes extreme. The coordination for the latter is usually called "stigmergy". Stigmergy is a generic name for mechanisms that provide spontaneous, indirect coordination between agents, where the influence in the environment left by the behavior of one agent stimulates the performance of a subsequent action of this agent or a different agent [7]. Since direct coordination is unnecessary in stigmergy, this mechanism can work in situations with massive numbers of agents. However, there is no guarantee that optimal coordination can be achieved. Therefore, how to create optimal coordination using stigmergy is an ambitious topic for research.

Moreover, traffic-control systems essentially have an interesting feature for the system architecture of MASs. Each agent in a MAS usually behaves to achieve a MAS goal regardless of its local or global views, and no agent behaves selfishly for its own gain. Of course, the goal of agents in a market-based environment is their own gain and basically they do behave selfishly. However, in the MAS for traffic-control systems, two competitive goals need to be achieved: the "goal of each agent" and the "goal of the MAS".

In the MAS for a traffic-control system each agent, which controls each car1, wants to behave selfishly to achieve its goal, e.g., optimal-route navigation by considering the shortest route and the avoidance of congestion. Therefore, each agent in the ITS is in a competitive situation similar to the conventional game environment in a MAS. However, the goal for the MAS itself is stability and optimizing the traffic-control system. That is, eliminating traffic congestion and minimizing the average travel time of all cars to attain a smooth traffic flow. To achieve these goals, it may be necessary to restrict the behavior of each agent. Consequently, the goal of each agent and the goal of the MAS have a competitive relation, and ITS is a very interesting application for the MAS.

In the VICS and PCS currently operated in Japan, congestion information and congestion-forecasting information are updated every 5 min. In other words, VISC or PCS cannot forecast less than five minutes ahead. Since traffic data from sensors and cars are collected at the central server and calculations are done by using all the data, this process

<sup>1</sup> A system that interacts with a human driver to lead him to a destination as in a car-navigation system.

needs a certain amount of time. Therefore, we propose a new congestion-forecasting system that can react to dynamically changing traffic conditions based on a coordination mechanism using the pheromone-communication model. Its main feature is to be able to forecast short-term congestion one or two minutes ahead. There have been many studies on ITS [8][9][10], but there have been few on the forecasting of short-term congestion.

Section 2 discusses traffic-congestion information for car agents and proposes a method of forecasting congestion that uses a multi-agent coordination mechanism for road agents set up at intersections based on a pheromone-communications model, which adaptively responds to increasing amounts of congestion. Section 3 discusses our tests to verify the basic effectiveness of this method. Finally, we conclude this paper in Section 4.

## **2. Method of forecasting congestion based on pheromonecommunications model**

2 Steroids

in these systems.

PCS also provides information on car movements and that on forecasting traffic congestion to individual drivers the same as VICS does. Different from VICS, this system collects traffic information from all cars, which are considered to be movable sensor units. Each car has a telecommunication unit and transmits several kinds of information such as position, velocity, and the status of the car to the central server. Then, the calculated car-movement and traffic-congestion-forecasting information are analyzed and the information from the center

Though these two systems are currently operated in Japan, the system structure of both systems is top down and centralized, so the reaction to dynamic changes in traffic congestion and occurrence of accidents is usually delayed and serious problems can occur when the central server is down. In other words, there is a lack of real-time features and of robustness

On the other hand, traffic-control systems like ITS and PCS essentially have interesting features for the coordination mechanisms of multi-agent systems (MASs). The coordination mechanisms of MAS can generally be divided into two types: direct and indirect. In the former, precise coordination can be achieved, but when the number of agents becomes excessive the load of coordination becomes extreme. The coordination for the latter is usually called "stigmergy". Stigmergy is a generic name for mechanisms that provide spontaneous, indirect coordination between agents, where the influence in the environment left by the behavior of one agent stimulates the performance of a subsequent action of this agent or a different agent [7]. Since direct coordination is unnecessary in stigmergy, this mechanism can work in situations with massive numbers of agents. However, there is no guarantee that optimal coordination can be achieved. Therefore, how to create optimal coordination using

Moreover, traffic-control systems essentially have an interesting feature for the system architecture of MASs. Each agent in a MAS usually behaves to achieve a MAS goal regardless of its local or global views, and no agent behaves selfishly for its own gain. Of course, the goal of agents in a market-based environment is their own gain and basically they do behave selfishly. However, in the MAS for traffic-control systems, two competitive goals need to be

In the MAS for a traffic-control system each agent, which controls each car1, wants to behave selfishly to achieve its goal, e.g., optimal-route navigation by considering the shortest route and the avoidance of congestion. Therefore, each agent in the ITS is in a competitive situation similar to the conventional game environment in a MAS. However, the goal for the MAS itself is stability and optimizing the traffic-control system. That is, eliminating traffic congestion and minimizing the average travel time of all cars to attain a smooth traffic flow. To achieve these goals, it may be necessary to restrict the behavior of each agent. Consequently, the goal of each agent and the goal of the MAS have a competitive relation, and ITS is a very

In the VICS and PCS currently operated in Japan, congestion information and congestion-forecasting information are updated every 5 min. In other words, VISC or PCS cannot forecast less than five minutes ahead. Since traffic data from sensors and cars are collected at the central server and calculations are done by using all the data, this process

<sup>1</sup> A system that interacts with a human driver to lead him to a destination as in a car-navigation system.

is displayed on equipment, such as car-navigation systems.

stigmergy is an ambitious topic for research.

interesting application for the MAS.

achieved: the "goal of each agent" and the "goal of the MAS".

Congestion-forecasting technology is one of the main elements of ITS. Up to now, several methods have been proposed, two of which are classified below.


Although it can effectively make forecasts under regular-congestion conditions that have originated from car and road situations, a large amount of past data is necessary for analysis. Moreover, it has weaknesses in forecasting under irregular-congestion conditions, such as those experienced during the Golden-week holidays in Japan and the Christmas-holiday season in the U.S.

VICS and PCS essentially belong to the second classification, and is excellent at short-term forecasting of congestion. Yet, in the current VICS and PCS that is operating in Japan, data from each car is collected at the central server and all calculations are done there. Consequently, it is difficult to supply real-time information due to bottlenecks and time lags in communicating information and the centralized calculations.

For solving these problems, a short-term system of forecasting congestion based on distributed processing is adequate. This paper focuses on "roads", and we propose a MAS consisting of many road agents. These road agents are set up at every intersection, and they coordinate locally with one another to forecast congestion. In this research, we adopted the pheromone-communications model as the mechanism for coordination.

Pheromone communications are based on the behavior of social insects like ants and bees and are applied as a model that is used to adaptively respond to dynamic changes in the environment in various applications [12] (see Fig. 2). Our method is an effective way of forecasting congestion in which each road agent generates pheromone information on its own road unit and exchanges this information with its neighboring agents. In a related study, Ando et al. investigated the forecasting of congestion in a local area a short time after pheromones had evaporated and diffused [11]. However, all drivers need to have the same probe-car system installed in their vehicles. Even though there has been some discussion on

When one ant finds an advantageous path from the colony to food, others are more likely to follow that path, and positive feedback eventually leads all the ants to follow a single path.

,

Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model

http://dx.doi.org/10.5772/52563

167

• A road unit is a section between two connected intersections. Each road unit consists of

• The number of cars going through an intersection is counted by a sensor installed at each intersection, and this number is sent to each road agent installed on roadside server

• The road agent installed in each roadside server computer calculates and forecasts the

A road unit on which a car is currently traveling is called "upstream", and a road unit that will be reached in the future is called "downstream". We focused on two important car-flow dynamics to investigate traffic congestion (see Fig. 4). The first was the flow in traffic density, which spreads from upstream to downstream, corresponding to the movement of cars. The second was the flow in traffic congestion, which spreads from downstream to upstream. At this point, the traffic congestion is defined as follows: a certain road unit becomes bottle-necked blocking the flow of cars. This blocking generates a queue of cars

Therefore, central servers and probe-car systems are not necessary with our method.

,

First, we will define the road environment as follows (see Fig. 3):

several lanes, usually in both directions, with no branching.

**Figure 5.** Flow chart for forecasting congestion

computers at regular intervals.

traffic congestion.

from downstream to upstream.

**2.1. Congestion-forecasting algorithm**

**Figure 2.** Ant-colony optimization

information being shared, individual automobile manufacturers are currently developing their own probe-car systems and consequently the rate of diffusion of these probe-car systems is quite low. Therefore, we decided to develop a more realistic and universal system by focusing on the road and not the cars.

**Figure 4.** Two important flows in congestion dynamics

<sup>166</sup> Ant Colony Optimization - Techniques and Applications Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model 5 Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model http://dx.doi.org/10.5772/52563 167

**Figure 5.** Flow chart for forecasting congestion

4 Steroids

When one ant finds an advantageous path from the colony to food, others are more likely to follow that path, and positive feedback eventually leads all the

information being shared, individual automobile manufacturers are currently developing their own probe-car systems and consequently the rate of diffusion of these probe-car systems is quite low. Therefore, we decided to develop a more realistic and universal system

ants to follow a single path.

by focusing on the road and not the cars.

**Figure 2.** Ant-colony optimization

**Figure 3.** Structure of road environment

**Figure 4.** Two important flows in congestion dynamics

## **2.1. Congestion-forecasting algorithm**

First, we will define the road environment as follows (see Fig. 3):


Therefore, central servers and probe-car systems are not necessary with our method.

A road unit on which a car is currently traveling is called "upstream", and a road unit that will be reached in the future is called "downstream". We focused on two important car-flow dynamics to investigate traffic congestion (see Fig. 4). The first was the flow in traffic density, which spreads from upstream to downstream, corresponding to the movement of cars. The second was the flow in traffic congestion, which spreads from downstream to upstream. At this point, the traffic congestion is defined as follows: a certain road unit becomes bottle-necked blocking the flow of cars. This blocking generates a queue of cars from downstream to upstream.

**Figure 6.** Calculation of current traffic situation

In this paper, we formulate the flow of traffic density using "traffic-density pheromones ∆*τ*" and formulate the growth of the queue using "congestion-diffusion pheromones *q*". To make forecasts more accurate, we introduce the "evaporation rate *e*", which indicates the change in congestion density from generation to dissolution.

Each road agent in our algorithm forecasts traffic congestion at one minute intervals, as shown in Fig. 5, where *τ*(*p*, *t*, *x*) is the forecasted traffic density of a road unit *s* at time *t* and ∆*τ*(*p*, *t*, *x*) is the forecasted transition in traffic congestion of road unit *s* at time *t*. Even though the calculation interval for forecasting can be shortened further, this increases the load of communication between the sensor and road agent. The one-minute intervals are much shorter than the five minutes of VICS.

Forecasting one minute ahead is calculated through coordination between each agent and their adjacent neighbouring agents. And forecasting two or more minutes ahead is calculated through coordination between each agent and more dispersed neighbouring agents.

## **2.2. Calculation of the current traffic situation**

The inflowing amount, *I*(*p*, *t*), and outflowing number, *O*(*p*, *t*), of cars at regular intervals *t* are measured with a sensor and are sent to road agents. *I*(*p*, *t*) indicates how many cars flowed into a road unit, *p*, and *O*(*p*, *t*) indicates how many flowed out of it. First, the road agent that receives this information calculates the traffic density as

$$N(p,t) = N(p, t-1) + I(p, t) - O(p, t) \tag{1}$$

$$d(p, t) = \frac{N(p, t) \times l\_{car}}{l\_p \times L\_p} \tag{2}$$

**Figure 7.** Traffic simulator

, *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>), and <sup>∆</sup>*τ*(*p*′

*<sup>τ</sup>*(*p*′

**2.3. Calculation of congestion forecasting pheromone**

from the neighbouring road unit. Then, *τ*(*p*, *t*, *x*) is calculated.

**(a) Calculation of traffic-density pheromones**

of traffic density due to the congestion.

Each road agent calculates the congestion forecasting pheromone, *τ*, which indicates the forecasted congestion density that will occur a few minutes ahead the current situation. *τ*(*p*, *t*, 0) = *d*(*p*, *t*) and ∆*τ*(*p*, *t*, 0) = *I*(*p*, *t*) − *O*(*p*, *t*) are the initial values for this calculation. As Fig. 5 shows, the traffic-density pheromones ∆*τ*(*p*, *t*, *x*), congestion-diffusion pheromones *q*(*p*, *t*, *x*), and evaporation rate *e*(*p*, *t*, *x*) are calculated using *τ*(*p*, *t*, *x* − 1), ∆*τ*(*p*, *t*, *x* − 1),

As previously mentioned, the traffic density spreads from upstream to downstream, corresponding to the movement of cars. What is important is how fast this flow is

Here, *sp* is the distance that a car moves during a certain time span and this is calculated from the maximum legal speed limit of road unit *p*. *bsp* is the proportion of time green lights are displayed in a signal cycle in the traveling direction of the car on this road. Moreover, *j f*(*p*, *t*, *x*) is a congestion factor that shows the decreasing ratio of the transmitting velocity

, *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>) and <sup>∆</sup>*τ*(*p*′

Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model

http://dx.doi.org/10.5772/52563

169

*S*(*p*, *t*, *x*) = *sp* × *bsp* × *j f*(*p*, *t*, *x*) (3)

, *t*, *x* − 1) are given

, *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>). At this point, *<sup>τ</sup>*(*p*′

transmitted, and we define the transmitting velocity of traffic density as *S*(*p*, *t*, *x*).

where *N*(*p*, *t*) is the number of cars, *d*(*p*, *t*) is the traffic density at intervals *t* of a road unit, *p*, *lcar* is the length of a car2, *lp* is the length of the road unit, *p*, and *Lp* is the number of lanes of *p* (see Fig. 6).

<sup>2</sup> More precisely, the length of a car + the distance between two cars.

**Figure 7.** Traffic simulator

6 Steroids

agents.

lanes of *p* (see Fig. 6).

**Figure 6.** Calculation of current traffic situation

in congestion density from generation to dissolution.

much shorter than the five minutes of VICS.

**2.2. Calculation of the current traffic situation**

<sup>2</sup> More precisely, the length of a car + the distance between two cars.

agent that receives this information calculates the traffic density as

In this paper, we formulate the flow of traffic density using "traffic-density pheromones ∆*τ*" and formulate the growth of the queue using "congestion-diffusion pheromones *q*". To make forecasts more accurate, we introduce the "evaporation rate *e*", which indicates the change

Each road agent in our algorithm forecasts traffic congestion at one minute intervals, as shown in Fig. 5, where *τ*(*p*, *t*, *x*) is the forecasted traffic density of a road unit *s* at time *t* and ∆*τ*(*p*, *t*, *x*) is the forecasted transition in traffic congestion of road unit *s* at time *t*. Even though the calculation interval for forecasting can be shortened further, this increases the load of communication between the sensor and road agent. The one-minute intervals are

Forecasting one minute ahead is calculated through coordination between each agent and their adjacent neighbouring agents. And forecasting two or more minutes ahead is calculated through coordination between each agent and more dispersed neighbouring

The inflowing amount, *I*(*p*, *t*), and outflowing number, *O*(*p*, *t*), of cars at regular intervals *t* are measured with a sensor and are sent to road agents. *I*(*p*, *t*) indicates how many cars flowed into a road unit, *p*, and *O*(*p*, *t*) indicates how many flowed out of it. First, the road

*<sup>d</sup>*(*p*, *<sup>t</sup>*) = *<sup>N</sup>*(*p*, *<sup>t</sup>*) <sup>×</sup> *lcar*

where *N*(*p*, *t*) is the number of cars, *d*(*p*, *t*) is the traffic density at intervals *t* of a road unit, *p*, *lcar* is the length of a car2, *lp* is the length of the road unit, *p*, and *Lp* is the number of

*lp* × *Lp*

*N*(*p*, *t*) = *N*(*p*, *t* − 1) + *I*(*p*, *t*) − *O*(*p*, *t*) (1)

(2)

## **2.3. Calculation of congestion forecasting pheromone**

Each road agent calculates the congestion forecasting pheromone, *τ*, which indicates the forecasted congestion density that will occur a few minutes ahead the current situation. *τ*(*p*, *t*, 0) = *d*(*p*, *t*) and ∆*τ*(*p*, *t*, 0) = *I*(*p*, *t*) − *O*(*p*, *t*) are the initial values for this calculation. As Fig. 5 shows, the traffic-density pheromones ∆*τ*(*p*, *t*, *x*), congestion-diffusion pheromones *q*(*p*, *t*, *x*), and evaporation rate *e*(*p*, *t*, *x*) are calculated using *τ*(*p*, *t*, *x* − 1), ∆*τ*(*p*, *t*, *x* − 1), *<sup>τ</sup>*(*p*′ , *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>), and <sup>∆</sup>*τ*(*p*′ , *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>). At this point, *<sup>τ</sup>*(*p*′ , *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>) and <sup>∆</sup>*τ*(*p*′ , *t*, *x* − 1) are given from the neighbouring road unit. Then, *τ*(*p*, *t*, *x*) is calculated.

## **(a) Calculation of traffic-density pheromones**

As previously mentioned, the traffic density spreads from upstream to downstream, corresponding to the movement of cars. What is important is how fast this flow is transmitted, and we define the transmitting velocity of traffic density as *S*(*p*, *t*, *x*).

$$S(p, t, \mathbf{x}) = s\_p \times bs\_p \times jf(p, t, \mathbf{x}) \tag{3}$$

Here, *sp* is the distance that a car moves during a certain time span and this is calculated from the maximum legal speed limit of road unit *p*. *bsp* is the proportion of time green lights are displayed in a signal cycle in the traveling direction of the car on this road. Moreover, *j f*(*p*, *t*, *x*) is a congestion factor that shows the decreasing ratio of the transmitting velocity of traffic density due to the congestion.

$$f(p, t, \mathbf{x}) = \begin{cases} 1.0 & (\tau(p, t, \mathbf{x} - 1) < a) \\ 1.0 - \tau(p, t, \mathbf{x} - 1) & (\tau(p, t, \mathbf{x} - 1) \ge a) \end{cases} \tag{4}$$

*α* is a threshold where the congestion factor demonstrates the effect, and *α* = 0.5 is used here. *S*(*p*, *t*, *x*) indicates the transmission distance of the traffic density in the one time span, so the ratio of *S*(*p*, *t*, *x*) and *lp* is important.

$$
\Delta\tau(p, t, \mathbf{x}) = \sum\_{p' \subset N\_b} f(p, p') \times \Delta\tau'(p', p, t, \mathbf{x}) \tag{5}
$$

**Figure 8.** Forecasting scenario in simulator

**Figure 9.** Comparison of congestion forecasting due to changes in traffic density (forecasting 1 min ahead).

decentralization, *vp*, of this change is calculated using the data from a previous day.

*e*(*p*, *t*, *x*) =

 

∆*τ*(*p*, *t*, 0) = *I*(*p*, *t*) − *O*(*p*, *t*), ∆*τ* becoming larger than normal means that traffic congestion will occur. However, ∆*τ* becoming smaller than normal means that traffic congestion is "evaporating". To determine the amount of normal traffic change on the road unit *p*, the

*<sup>β</sup>*1&(*v*(*p*, *<sup>t</sup>*, 0) > *<sup>x</sup>* × *vp*)

*<sup>β</sup>*2&(*v*(*p*, *<sup>t</sup>*, 0) < −*<sup>x</sup>* × *vp*)

1.0&(−*x* × *vp* < *v*(*p*, *t*, 0) < *x* × *vp*)

Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model

http://dx.doi.org/10.5772/52563

171

(9)

$$\Delta \pi'(p, t, \mathbf{x}) = \begin{cases} \Delta \pi(p, t, \mathbf{x} - 1) \& (S(p, t, \mathbf{x}) > l\_p) \\ \frac{S(p, t, \mathbf{x})}{l\_p} \times \Delta \pi(p, t, \mathbf{x} - 1) \& (S(p, t, \mathbf{x}) \le l\_p) \end{cases} \tag{6}$$

where *Nb* indicates the set of upstream road units *<sup>p</sup>* and *<sup>f</sup>*(*p*, *<sup>p</sup>*′ ) is a parameter that changes based on the relation between *<sup>p</sup>* and *<sup>p</sup>*′ . In this study, *<sup>f</sup>*(*p*, *<sup>p</sup>*′ ) was 0.7 when the road unit, *<sup>p</sup>*′ → *<sup>p</sup>*, was straight, 0.2 when it turned left, and 0.1 when it turned right.

#### **(b) Calculation of congestion diffusion pheromones**

As previously mentioned, traffic congestion spreads from downstream to upstream. Therefore, the congestion diffusion pheromones are defined based on the difference between the congestion level of the current road unit and the congestion level of the next road unit.

$$q(p, t, \mathbf{x}) = \sum\_{p'' \subset N\_f} \mathbf{g}(p, p'') \times q'(p'', p, t, \mathbf{x}) \tag{7}$$

$$q'(p'',p,t,\mathbf{x}) = \{\mathbf{\tau}(p'',t,\mathbf{x}-1) - \mathbf{\tau}(p,t,\mathbf{x}-1)\}\tag{8}$$

where *Nf* indicates the set of downstream road unit *p*, and *g*(*p*, *p*") is a parameter that changes based on the relation between *<sup>p</sup>* and *<sup>p</sup>*′′, which is the same as *<sup>f</sup>*(*p*, *<sup>p</sup>*′ ).


**Table 1.** Comparison of correlation coefficient between forecast and actual values due to changes in traffic density (forecasting 1 min ahead).

#### **(c) Calculation of evaporation rate**

As previously mentioned, the evaporation rate indicates the change in congestion density due to its generation and dissolution. That is, by referring to the degree of traffic change, <sup>170</sup> Ant Colony Optimization - Techniques and Applications Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model 9 Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model http://dx.doi.org/10.5772/52563 171

**Figure 8.** Forecasting scenario in simulator

8 Steroids

*j f*(*p*, *t*, *x*) =

<sup>∆</sup>*τ*(*p*, *<sup>t</sup>*, *<sup>x</sup>*) = ∑

where *Nb* indicates the set of upstream road units *<sup>p</sup>* and *<sup>f</sup>*(*p*, *<sup>p</sup>*′

*<sup>q</sup>*(*p*, *<sup>t</sup>*, *<sup>x</sup>*) = ∑

*p*"⊂*Nf*

changes based on the relation between *<sup>p</sup>* and *<sup>p</sup>*′′, which is the same as *<sup>f</sup>*(*p*, *<sup>p</sup>*′

*S*(*p*,*t*,*x*)

*<sup>p</sup>*′ → *<sup>p</sup>*, was straight, 0.2 when it turned left, and 0.1 when it turned right.

*p*′⊂*Nb*

so the ratio of *S*(*p*, *t*, *x*) and *lp* is important.

∆*τ*′

based on the relation between *<sup>p</sup>* and *<sup>p</sup>*′

(*p*, *t*, *x*) =

**(b) Calculation of congestion diffusion pheromones**

*q*′

(forecasting 1 min ahead).

**(c) Calculation of evaporation rate**

1.0 (*τ*(*p*, *t*, *x* − 1) < *α*)

) <sup>×</sup> <sup>∆</sup>*τ*′

(*p*′

*α* is a threshold where the congestion factor demonstrates the effect, and *α* = 0.5 is used here. *S*(*p*, *t*, *x*) indicates the transmission distance of the traffic density in the one time span,

*<sup>f</sup>*(*p*, *<sup>p</sup>*′

∆*τ*(*p*, *t*, *x* − 1)&(*S*(*p*, *t*, *x*) > *lp*)

As previously mentioned, traffic congestion spreads from downstream to upstream. Therefore, the congestion diffusion pheromones are defined based on the difference between the congestion level of the current road unit and the congestion level of the next road unit.

where *Nf* indicates the set of downstream road unit *p*, and *g*(*p*, *p*") is a parameter that

Forecasting with pheromone method 0.98 0.94 0.91 Conventional forecasting 0.95 0.88 0.79

As previously mentioned, the evaporation rate indicates the change in congestion density due to its generation and dissolution. That is, by referring to the degree of traffic change,

**Table 1.** Comparison of correlation coefficient between forecast and actual values due to changes in traffic density

*<sup>g</sup>*(*p*, *<sup>p</sup>*") <sup>×</sup> *<sup>q</sup>*′

(*p*", *p*, *t*, *x*) = {*τ*(*p*", *t*, *x* − 1) − *τ*(*p*, *t*, *x* − 1)} (8)

**1 min ahead 3 min ahead 5 min ahead**

. In this study, *<sup>f</sup>*(*p*, *<sup>p</sup>*′

1.0 <sup>−</sup> *<sup>τ</sup>*(*p*, *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>) (*τ*(*p*, *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>) <sup>≥</sup> *<sup>α</sup>*) (4)

*lp* <sup>×</sup> <sup>∆</sup>*τ*(*p*, *<sup>t</sup>*, *<sup>x</sup>* <sup>−</sup> <sup>1</sup>)&(*S*(*p*, *<sup>t</sup>*, *<sup>x</sup>*) <sup>≤</sup> *lp*) (6)

, *p*, *t*, *x*) (5)

) is a parameter that changes

) was 0.7 when the road unit,

(*p*", *p*, *t*, *x*) (7)

).

**Figure 9.** Comparison of congestion forecasting due to changes in traffic density (forecasting 1 min ahead).

∆*τ*(*p*, *t*, 0) = *I*(*p*, *t*) − *O*(*p*, *t*), ∆*τ* becoming larger than normal means that traffic congestion will occur. However, ∆*τ* becoming smaller than normal means that traffic congestion is "evaporating". To determine the amount of normal traffic change on the road unit *p*, the decentralization, *vp*, of this change is calculated using the data from a previous day.

$$e(p, t, \mathbf{x}) = \begin{cases} \beta\_1 \pounds(v(p, t, 0) > \mathbf{x} \times \mathbf{v}\_p) \\ 1.0 \pounds(-\mathbf{x} \times \mathbf{v}\_p < v(p, t, 0) < \mathbf{x} \times \mathbf{v}\_p) \\ \beta\_2 \pounds(v(p, t, 0) < -\mathbf{x} \times \mathbf{v}\_p) \end{cases} \tag{9}$$

The above expression shows that when the difference between the observed amount of traffic and the amount of traffic in normal conditions increases, the degree of congestion generation and congestion evaporation becomes big. *β*<sup>1</sup> and *β*<sup>2</sup> are parameters that indicate the degree of evaporation and in this study, *β*<sup>1</sup> is 1.1 and *β*<sup>2</sup> is 0.9.

## **(d) Calculating congestion-forecasting pheromones**

The forecasting of congestion pheromones after *x* minutes is calculated from the above value as follows.

$$
\tau(p, t, \mathbf{x}) = e(p, t, \mathbf{x}) \times \tau(p, t, \mathbf{x} - 1) + \Delta \tau(p, t, \mathbf{x}) + q(p, t, \mathbf{x}) \tag{10}
$$

**Figure 10.** Comparison of congestion forecasts due to sudden accidents (forecasting 1 min ahead)

road unit that will continue for longer than 5 min.

*2.4.1. Congestion due to changes in traffic density*

*2.4.2. Congestion due to sudden accidents*

congestion more accurately than the conventional approach.

occur after this. On the other hand, we can also see that one agent forecasts congestion of its

Traffic congestion is usually generated when more than the acceptable number of cars moves into a road unit. We carried out the simulation for about 2 hours and generated and evaporated congestion several times by changing the traffic density. We then compared our proposed method with the conventional approach by forecasting 1, 3, and 5 min ahead.

As a result, our proposed approach had a higher accuracy than the conventional method (Table 1). Fig. 9 shows the change in the actual traffic-congestion level (blue line) and the forecast congestion level 1 min ahead by using the conventional (yellow line) and our approach (red line). The change in the red line is similar to that in the blue line. The change in the yellow line, on the other hand, is delayed. Therefore, our proposed method can forecast

Next, we evaluated how accurately congestion could be forecast when sudden accidents occurred. This type of congestion does not happen based on changes in traffic density, but it occurs due to the decreased capacity of the roads to accommodate cars traveling along them. Since this decrease in capacity happens suddenly, the speed at which congestion is diffused is very rapid. In our simulation, we compared the effectiveness of our proposed method with that of the conventional approach by quickly changing the traffic density of a

Forecasting with pheromone method 0.98 0.94 0.86 Conventional forecasting 0.86 0.66 0.45

**Table 2.** Comparison of correlation coefficient between forecast and actual values due to sudden accidents

**1 min ahead 3 min ahead 5 min ahead**

Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model

http://dx.doi.org/10.5772/52563

173

Each road agent forecasts short-term traffic congestion by sequentially and repeatedly calculating (a) to (d) from the above.

## **2.4. Simulations**

To experimentally verify the basic effectiveness of our proposed forecasting model, we implemented a simple simulation environment and compared the accuracy of forecasting a few minutes ahead (i.e., one, three, and five minutes) with the proposed and a conventional method. We especially verified the effectiveness of our methodology in two respects, i.e.,


The correlation coefficient of the actual measurements and the forecasting values was used for the evaluation, and the simulation environment shown in Fig. 7 was used for the experiment. This simple simulator had a 5 x 5 lattice structure with single-lane roads. There was one traffic signal at each intersection. The length of one road unit, i.e., the distance between two consecutive intersections, was 400 m.

In this experiment, we set 1*time* − *step* to 1*sec* and 1*time* − *span* to 60*time* − *steps*. Each road agent calculated the traffic density on its own road unit every minute, and forecasted until 5 minutes ahead. To evaluate the effectiveness of the proposed method, we used a conventional method of short-term forecasting based on a statistical approach [8] and made forecasts 1, 3, and 5 min ahead.

This conventional forecasting approach was based on the assumption that the current congestion situation would generally continue for a few minutes. The current VICS and PCS update their congestion information every five minutes, so if we used this conventional method to forecast five minutes ahead, it could basically be thought of as using the same approach as VICS and PCS.

Fig. 8 is an expansion of part of the simulator used in executing the forecasts. We can see that one road agent forecasts congestion of its road unit that will not occur within 5 min but <sup>172</sup> Ant Colony Optimization - Techniques and Applications Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model 11 Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model http://dx.doi.org/10.5772/52563 173

**Figure 10.** Comparison of congestion forecasts due to sudden accidents (forecasting 1 min ahead)

occur after this. On the other hand, we can also see that one agent forecasts congestion of its road unit that will continue for longer than 5 min.


**Table 2.** Comparison of correlation coefficient between forecast and actual values due to sudden accidents

## *2.4.1. Congestion due to changes in traffic density*

10 Steroids

as follows.

**2.4. Simulations**

and 5 min ahead.

approach as VICS and PCS.

The above expression shows that when the difference between the observed amount of traffic and the amount of traffic in normal conditions increases, the degree of congestion generation and congestion evaporation becomes big. *β*<sup>1</sup> and *β*<sup>2</sup> are parameters that indicate

The forecasting of congestion pheromones after *x* minutes is calculated from the above value

Each road agent forecasts short-term traffic congestion by sequentially and repeatedly

To experimentally verify the basic effectiveness of our proposed forecasting model, we implemented a simple simulation environment and compared the accuracy of forecasting a few minutes ahead (i.e., one, three, and five minutes) with the proposed and a conventional method. We especially verified the effectiveness of our methodology in two respects, i.e.,

The correlation coefficient of the actual measurements and the forecasting values was used for the evaluation, and the simulation environment shown in Fig. 7 was used for the experiment. This simple simulator had a 5 x 5 lattice structure with single-lane roads. There was one traffic signal at each intersection. The length of one road unit, i.e., the distance

In this experiment, we set 1*time* − *step* to 1*sec* and 1*time* − *span* to 60*time* − *steps*. Each road agent calculated the traffic density on its own road unit every minute, and forecasted until 5 minutes ahead. To evaluate the effectiveness of the proposed method, we used a conventional method of short-term forecasting based on a statistical approach [8] and made forecasts 1, 3,

This conventional forecasting approach was based on the assumption that the current congestion situation would generally continue for a few minutes. The current VICS and PCS update their congestion information every five minutes, so if we used this conventional method to forecast five minutes ahead, it could basically be thought of as using the same

Fig. 8 is an expansion of part of the simulator used in executing the forecasts. We can see that one road agent forecasts congestion of its road unit that will not occur within 5 min but

*τ*(*p*, *t*, *x*) = *e*(*p*, *t*, *x*) × *τ*(*p*, *t*, *x* − 1) + ∆*τ*(*p*, *t*, *x*) + *q*(*p*, *t*, *x*) (10)

the degree of evaporation and in this study, *β*<sup>1</sup> is 1.1 and *β*<sup>2</sup> is 0.9.

**(d) Calculating congestion-forecasting pheromones**

calculating (a) to (d) from the above.

1. The forecasting accuracy of generation/dis-

2. The forecasting accuracy of generation/dissolution of congestion due to sudden accidents.

between two consecutive intersections, was 400 m.

solution of congestion due to changes in traffic density and

Traffic congestion is usually generated when more than the acceptable number of cars moves into a road unit. We carried out the simulation for about 2 hours and generated and evaporated congestion several times by changing the traffic density. We then compared our proposed method with the conventional approach by forecasting 1, 3, and 5 min ahead.

As a result, our proposed approach had a higher accuracy than the conventional method (Table 1). Fig. 9 shows the change in the actual traffic-congestion level (blue line) and the forecast congestion level 1 min ahead by using the conventional (yellow line) and our approach (red line). The change in the red line is similar to that in the blue line. The change in the yellow line, on the other hand, is delayed. Therefore, our proposed method can forecast congestion more accurately than the conventional approach.

## *2.4.2. Congestion due to sudden accidents*

Next, we evaluated how accurately congestion could be forecast when sudden accidents occurred. This type of congestion does not happen based on changes in traffic density, but it occurs due to the decreased capacity of the roads to accommodate cars traveling along them. Since this decrease in capacity happens suddenly, the speed at which congestion is diffused is very rapid. In our simulation, we compared the effectiveness of our proposed method with that of the conventional approach by quickly changing the traffic density of a certain road unit. As Table 2 and Fig. 10 show, our proposed method is more accurate than the conventional scheme. Forecasting accuracy particularly worsened with the conventional method during long-term forecasts. However, we were able to maintain accurate forecasts with our method.

[3] X.-H. Yu and W. W. Recker: Stochastic adaptive control model for traffic signal systems,

Traffic-Congestion Forecasting Algorithm Based on Pheromone Communication Model

http://dx.doi.org/10.5772/52563

175

[4] K. Tufte, J. Li, D. Maier, V. Papadimos, R. L. Bertini, and J. Rucker: Travel Time

[5] J.-D. Schmocker, S. Ahuja, and M. G. H. Bell: Multi-objective signal control of urban junctions - Framework and a London case study, Transportation Research Part C 16, pp.

[7] "Definitions of stigmergy." From a special Issue of Artificial Life on Stigmergy. Volume

[8] P. G. Balaji, G. Sachdeva, D. Srinivasan, Multi-agent System based Urban Traffic Management, IEEE Congress on Evolutionary Computation 2007, pp. 1740–1747, 2007.

[9] Ana L. C. Bazzan, Opportunities for multiagent systems and multiagent reinforcement learning in traffic control, Autonomous Agents and Multi-Agent Systems, Vol. 18, No.

[10] J. J. Sánchez-Medina and M. J. Gálan-Moreno and E. Rubio-Royo", Traffic Signal Optimization in La Almozara District in Saragossa Under Congestion Conditions, Using Genetic Algorithms, Traffic Microsimulation, and Cluster Computing, IEEE Transaction

[11] Y. Ando, Y. Fukazawa, O. Masutani, H. Iwasaki, and S. Honiden: Performance of Pheromone Model for Predicting Traffic Congestion, The Fifth International Joint Conference on Autonomous Agents and Multiagent Systems (AAMAS-06), 2006.

[12] M. Dorigo and G. di Caro: The ant colony optimization meta-heuristic, New ideas in

[13] E. Chung: Classification of traffic pattern. Proceedings of the 11th World Congress on

[14] M. Wiering: Multi-Agent Reinforcement Learning for Traffic Light Control, Machine Learning, *Proceedings of the Seventeenth International Conference (ICML'2000)*, pp.

[15] D. Oliveira, A. L. C. Bazzan, B. C. Silva, E. W. Basso, L. Nunes, and R. J. F. Rossetti, E. C. Oliveira, R. Silva, and L. C. Lamb: Reinforcement learning based control of traffic lights in non-stationary environments: A case study in a microscopic simulator, *Proc. of*

[16] D. Oliveira, P. Ferreira, and A. L. C. Bazzan: Reducing traffic jams with a swarm-based approach for selection of signal plans, *Proc. ANTS 2004*, Vol. 3172 of LNCS, pp. 416-417,

on Intelligent Transportation Systems, Vol. 11, No. 1. pp. 132–141, 2010.

Estimation Using NiagaraST and latte, SIGMOD'07, pp. 1091-1093, 2007.

Transportation Research Part C 14, pp. 263-282, Elsevier, 2006.

454-470, 2008.

5, Issue 2 / Spring 1999.

3, pp. 313–341, 2009.

ITS, 2003.

1151-1158, 2000.

*EUMAS06*, pp. 31-42, 2006.

Berlin, Germany, 2004.

[17] http://uscan.osoite.jp/

[6] http://www.vics.or.jp/english/index.html

optimization, McGrawHill, pp. 11-32, 1999.

## **3. Conclusion**

We proposed a method of forecasting congestion using a multi-agent coordination mechanism. A road agent installed at each intersection coordinates with its neighboring road agents based on the pheromone-communications model to adaptively respond to dynamically arising congestion and forecasts congestion a few minutes ahead. Here, we tested and verified the basic effectiveness of this method using simple simulation.

It is unnecessary in our approach to utilize a sufficient number of cars with the same probe system [8], or to upgrade the central server. At the very least, it needs to have simple sensors installed to count the number of cars moving through intersections, and small computers for road agents at these intersections. However, we have assumed that various kinds of computers, servers, and sensors will be installed in various locations to gather large amounts of information from the real world in about 5-10 years as part of urban scanning. Actually, small-scale real-world experiments are now being conducted in several locations throughout Japan [17]. These are based on the development of ubiquitous-information-communication technologies such as sensor-networks and wireless communication devices. In such situations, our method is expected to be quite practical. This evaluation was only done through simulation, and the road map used had a simple lattice structure. However, as we have already obtained detailed road data and VICS/PCS data throughout the entire country of Japan, we can shortly begin to evaluate our method using these real-world data.

As for traffic light control, all traffic lights need to react to dynamic changes in traffic in real time [14][15][16] in traffic-light-control systems, which are also the primary systems for controlling traffic. However, the current system cannot respond to dynamic changes in road conditions in real time even though some automatic control occurs according to the traffic flow. We also plan to develop a new traffic-light-control algorithm based on a MAS.

## **Author details**

Satoshi Kurihara

Osaka University, Japan

## **References**


12 Steroids

with our method.

**3. Conclusion**

using these real-world data.

**Author details**

Satoshi Kurihara

**References**

Osaka University, Japan

Multiagent Systems (AAMAS-06), 2006.

Prediction, KDD'04, pp. 22-25, 2004.

certain road unit. As Table 2 and Fig. 10 show, our proposed method is more accurate than the conventional scheme. Forecasting accuracy particularly worsened with the conventional method during long-term forecasts. However, we were able to maintain accurate forecasts

We proposed a method of forecasting congestion using a multi-agent coordination mechanism. A road agent installed at each intersection coordinates with its neighboring road agents based on the pheromone-communications model to adaptively respond to dynamically arising congestion and forecasts congestion a few minutes ahead. Here, we

It is unnecessary in our approach to utilize a sufficient number of cars with the same probe system [8], or to upgrade the central server. At the very least, it needs to have simple sensors installed to count the number of cars moving through intersections, and small computers for road agents at these intersections. However, we have assumed that various kinds of computers, servers, and sensors will be installed in various locations to gather large amounts of information from the real world in about 5-10 years as part of urban scanning. Actually, small-scale real-world experiments are now being conducted in several locations throughout Japan [17]. These are based on the development of ubiquitous-information-communication technologies such as sensor-networks and wireless communication devices. In such situations, our method is expected to be quite practical. This evaluation was only done through simulation, and the road map used had a simple lattice structure. However, as we have already obtained detailed road data and VICS/PCS data throughout the entire country of Japan, we can shortly begin to evaluate our method

As for traffic light control, all traffic lights need to react to dynamic changes in traffic in real time [14][15][16] in traffic-light-control systems, which are also the primary systems for controlling traffic. However, the current system cannot respond to dynamic changes in road conditions in real time even though some automatic control occurs according to the traffic

[1] V. R. Tomás and L. A. Garcia: A Cooperative Multiagent System for Traffic Management and Control, The Fourth International Joint Conference on Autonomous Agents and

[2] T. Nakata and J. Takeuchi: Mining Traffic Data from Probe-Car System for Travel Time

flow. We also plan to develop a new traffic-light-control algorithm based on a MAS.

tested and verified the basic effectiveness of this method using simple simulation.


**Chapter 8**

**Ant Colony Algorithm with**

E. H. Hay and K. Bertrand

http://dx.doi.org/10.5772/52051

of the features performance.

**1. Introduction**

**Applications in the Field of Genomics**

R. Rekaya, K. Robbins, M. Spangler, S. Smith,

Additional information is available at the end of the chapter

Ant colony algorithms (ACA) were first proposed by Dorigo *et al*. (1999) to solve difficult optimization problems, such as the traveling salesman, and have since been extended to solve many discrete optimization problems. As the name would imply, ACA are derived from the process by which ant colonies find the shortest route to a food source. Real ant col‐ onies communicate through the use of chemicals called pheromones which are deposited along the path an ant travels. Ants that choose a shorter path will transverse the distance at a faster rate, thus depositing more pheromone. Subsequent ants will then choose the path with more pheromone creating a positive feedback system. Artificial ants work as parallel units that communicate through a cumulative distribution function (CDF) that is updated by weights, determined by the "distance" traveled on a selected "path", which are analo‐ gous to the pheromones deposited by real ants (Dorigo *et al*. 1999, Ressom *et al.* 2006). As the CDF is updated, "paths" that perform better will be sampled at higher likelihoods by subse‐ quent artificial ants which, in turn, deposit more "pheromone", thus leading to a positive feedback system similar to the method of communication observed in real ant colonies. In the specific application of feature selection, the "path" chosen by an artificial ant is a subset of features selected from a larger sample space, and the "distance" traveled is some measure

The idea of selecting a sub-set of features capable of best classifying a group of samples can be, and has been, viewed as an optimization problem. The genetic algorithm (GA), simulated annealing (SA), and other optimization and machine learning algorithms have been applied to the problem of feature selection (Lin et al., 2006; Ooi and Tan, 2003; Peng et al., 2003; Albrecht et al., 2003). Though these methods are powerful, when deal‐

> © 2013 Rekaya et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

© 2013 Rekaya et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

distribution, and reproduction in any medium, provided the original work is properly cited.

**Chapter 8**

## **Ant Colony Algorithm with Applications in the Field of Genomics**

R. Rekaya, K. Robbins, M. Spangler, S. Smith, E. H. Hay and K. Bertrand

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/52051

## **1. Introduction**

Ant colony algorithms (ACA) were first proposed by Dorigo *et al*. (1999) to solve difficult optimization problems, such as the traveling salesman, and have since been extended to solve many discrete optimization problems. As the name would imply, ACA are derived from the process by which ant colonies find the shortest route to a food source. Real ant col‐ onies communicate through the use of chemicals called pheromones which are deposited along the path an ant travels. Ants that choose a shorter path will transverse the distance at a faster rate, thus depositing more pheromone. Subsequent ants will then choose the path with more pheromone creating a positive feedback system. Artificial ants work as parallel units that communicate through a cumulative distribution function (CDF) that is updated by weights, determined by the "distance" traveled on a selected "path", which are analo‐ gous to the pheromones deposited by real ants (Dorigo *et al*. 1999, Ressom *et al.* 2006). As the CDF is updated, "paths" that perform better will be sampled at higher likelihoods by subse‐ quent artificial ants which, in turn, deposit more "pheromone", thus leading to a positive feedback system similar to the method of communication observed in real ant colonies. In the specific application of feature selection, the "path" chosen by an artificial ant is a subset of features selected from a larger sample space, and the "distance" traveled is some measure of the features performance.

The idea of selecting a sub-set of features capable of best classifying a group of samples can be, and has been, viewed as an optimization problem. The genetic algorithm (GA), simulated annealing (SA), and other optimization and machine learning algorithms have been applied to the problem of feature selection (Lin et al., 2006; Ooi and Tan, 2003; Peng et al., 2003; Albrecht et al., 2003). Though these methods are powerful, when deal‐

© 2013 Rekaya et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Rekaya et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

ing with thousands of features across multiple classes, the computational cost of these methods can be prohibitive. Previous results obtained with these methods when dealing with large numbers of features, utilized filters to reduce the dimension of the datasets prior to implementation (Lin et al., 2006; Peng et al., 2006), or have produced relatively low prediction accuracies (Hong and Cho, 2006). For ACA, the communication of the ants through a common memory has a synergistic effect that, when coupled with more efficient searching of the sample space though the use of prior information, results in op‐ timal solutions being reached in far fewer iterations than required for GA or SA (Dorigo and Gambardella, 1997). The algorithm also lends itself to parallelization, with ants be‐ ing run on multiple processors, which can further reduce computation time, making its use more feasible with high dimension data sets.

## **2. General presentation of ant colony algorithm**

The ACA employs artificial ants that communicate through a probability density function (PDF) that is updated at-each iteration with weights or "pheromone levels", which are anal‐ ogous to the chemical pheromones used by real ants. The weights can be determined by the strength of the association between selected feature and the response of interest. Using the notation in [Dorigo and Gambardella, 1997; Ressom et al., 2006], the probability of sampling feature *m* at time *t* is defined as:

$$P\_m(t) = \frac{(\tau\_m(t))^\alpha \eta\_m{}^\beta}{\sum\_{m=1}^{\eta\_f^f} (\tau\_m(t))^\alpha \eta\_m{}^\beta} \tag{1}$$

Although the general idea of the ACA is simple and intuitive, its application to solve re‐ al world applications requires some good heuristics in defining the pheromone functions and their updating. In this chapter, we are presenting three applications of the ACA in the field of genetics and genomics based on previously published research by our group [Robbins et al., 2007, Robbins et al., 2008; Spangler et al., 2008; Rekaya and Robbins, 2009; Robbins et al., 2011]. Specific implementation details for each application are added

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

179

**2.1. Ant colony algorithm for feature selection in high dimension gene expression data for**

The idea of using gene expression data for diagnosis and personalized treatment presents a promising area of medicine and, as such, has been the focus of much research (Bagirov et al., 2003; Golub et al., 1999, Ramaswamy et al., 2001). Many algorithms have been developed to classify disease types based on the expression of selected genes, and significant gains have been made in the accuracy of disease classification (Antonov et al., 2004; Bagirov et al., 2003). In addition to the development of classification algorithms, many studies have shown that improved performance can be achieved when using a se‐ lected subset of features, as opposed to using all available data (Peng et al., 2003; Shen et al., 2006; Subramani et al., 2006). Increases in accuracy achieved through the selection of predictive features can complement and enhance the performance of classification al‐ gorithms, as well as improve the understanding of disease classes by identifying a small

In this section the ACA was implemented using the high-dimensional GCM data-set (Ram‐ aswamy et al., 2001), containing 16,063 genes and 14 tumor classes, with very limited prefiltering, and compared to several other rank based feature selection methods, as well as

*A.1 Latent variable model*: A Bayesian regression model was used to predict tumor type in the form of a probability *p*ic(yic=1), with yic = 1 indicating that sample i is from tumor class c. The regression on the vector of binary responses yc was done using a latent variable model (LVM), with *l*ic being an unobserved, continuous latent variable relating to binary response

where Xic corresponds to row i of the design matrix Xc for tumor class c. The link function of the expectation of the liability *Xicβc*with the binary response yic was constructed via a probit

previously published results to determine its efficacy as a feature selection method.

in the appropriate sections of the chapter.

set of biologically relevant features (Golub et al., 1999).

The liability *l*ic was modeled using a linear regression model as:

*ic*)= *Xicβ<sup>c</sup> eic* ~ *N* (0, 1)

model (West, 2003) yielding the following equations:

*pic*(*y*ic=1)=*Φ***(***Xicβc***)** and *pic*(*y*ic=1)=1−*Φ***(***Xicβc***)**

**disease classification**

yic such that:

1 if *l*

0 if *l*

*ic* = *Xicβ<sup>c</sup>* + *eic E*(*l*

*ic* ≥0

*ic* <0

*yic* ={

*l*

where *τm*(*t*) is the amount of pheromone for feature *m* at time *t*; *ηm*is some form of prior information on the expected performance of feature,*α* and *β* are parameters determining the weight given to pheromone deposited by ants and a priori information on the fea‐ tures, respectively.

Using the PDF as defined in equation (1), each of *j* artificial ants will select a subset *Sk* of *n* features from the sample space*S* containing all features. The pheromone level of each fea‐ ture *m* in *Sk* is then updated according to the performance of *Sk*as:

$$
\pi\_m(t+1) = (1-\rho)^\* \pi\_m(t) + \Delta \tau\_m(t) \tag{2}
$$

where *ρ* is a constant between 0 and 1 representing the rate at which the pheromone trail evaporates; *Δτm*(*t*)is the change in pheromone level for feature *m* based on the sum of accu‐ racy of all *Sk* containing SNP *m*, and is set to zero if feature *m* was not selected by any of the artificial ants.

Although the general idea of the ACA is simple and intuitive, its application to solve re‐ al world applications requires some good heuristics in defining the pheromone functions and their updating. In this chapter, we are presenting three applications of the ACA in the field of genetics and genomics based on previously published research by our group [Robbins et al., 2007, Robbins et al., 2008; Spangler et al., 2008; Rekaya and Robbins, 2009; Robbins et al., 2011]. Specific implementation details for each application are added in the appropriate sections of the chapter.

## **2.1. Ant colony algorithm for feature selection in high dimension gene expression data for disease classification**

The idea of using gene expression data for diagnosis and personalized treatment presents a promising area of medicine and, as such, has been the focus of much research (Bagirov et al., 2003; Golub et al., 1999, Ramaswamy et al., 2001). Many algorithms have been developed to classify disease types based on the expression of selected genes, and significant gains have been made in the accuracy of disease classification (Antonov et al., 2004; Bagirov et al., 2003). In addition to the development of classification algorithms, many studies have shown that improved performance can be achieved when using a se‐ lected subset of features, as opposed to using all available data (Peng et al., 2003; Shen et al., 2006; Subramani et al., 2006). Increases in accuracy achieved through the selection of predictive features can complement and enhance the performance of classification al‐ gorithms, as well as improve the understanding of disease classes by identifying a small set of biologically relevant features (Golub et al., 1999).

In this section the ACA was implemented using the high-dimensional GCM data-set (Ram‐ aswamy et al., 2001), containing 16,063 genes and 14 tumor classes, with very limited prefiltering, and compared to several other rank based feature selection methods, as well as previously published results to determine its efficacy as a feature selection method.

*A.1 Latent variable model*: A Bayesian regression model was used to predict tumor type in the form of a probability *p*ic(yic=1), with yic = 1 indicating that sample i is from tumor class c. The regression on the vector of binary responses yc was done using a latent variable model (LVM), with *l*ic being an unobserved, continuous latent variable relating to binary response yic such that:

$$y\_{ic} = \begin{cases} 1 & \text{if } l\_{ic} \ge 0 \\ 0 & \text{if } l\_{ic} < 0 \end{cases}$$

ing with thousands of features across multiple classes, the computational cost of these methods can be prohibitive. Previous results obtained with these methods when dealing with large numbers of features, utilized filters to reduce the dimension of the datasets prior to implementation (Lin et al., 2006; Peng et al., 2006), or have produced relatively low prediction accuracies (Hong and Cho, 2006). For ACA, the communication of the ants through a common memory has a synergistic effect that, when coupled with more efficient searching of the sample space though the use of prior information, results in op‐ timal solutions being reached in far fewer iterations than required for GA or SA (Dorigo and Gambardella, 1997). The algorithm also lends itself to parallelization, with ants be‐ ing run on multiple processors, which can further reduce computation time, making its

The ACA employs artificial ants that communicate through a probability density function (PDF) that is updated at-each iteration with weights or "pheromone levels", which are anal‐ ogous to the chemical pheromones used by real ants. The weights can be determined by the strength of the association between selected feature and the response of interest. Using the notation in [Dorigo and Gambardella, 1997; Ressom et al., 2006], the probability of sampling

1

t

( ( )) *m m*

=

where *τm*(*t*) is the amount of pheromone for feature *m* at time *t*; *ηm*is some form of prior information on the expected performance of feature,*α* and *β* are parameters determining the weight given to pheromone deposited by ants and a priori information on the fea‐

Using the PDF as defined in equation (1), each of *j* artificial ants will select a subset *Sk* of *n* features from the sample space*S* containing all features. The pheromone level of each fea‐

where *ρ* is a constant between 0 and 1 representing the rate at which the pheromone trail evaporates; *Δτm*(*t*)is the change in pheromone level for feature *m* based on the sum of accu‐ racy of all *Sk* containing SNP *m*, and is set to zero if feature *m* was not selected by any of the

( 1) (1 ) \* ( ) ( ) *<sup>m</sup> m m*

 rt

*m m m*

*t*

a b

 h

> t

*t tt* + = - +D (2)

å (1)

a b

 h

( ( )) ( )

t

*m nf*

ture *m* in *Sk* is then updated according to the performance of *Sk*as:

t

=

*<sup>t</sup> P t*

use more feasible with high dimension data sets.

178 Ant Colony Optimization - Techniques and Applications

feature *m* at time *t* is defined as:

tures, respectively.

artificial ants.

**2. General presentation of ant colony algorithm**

The liability *l*ic was modeled using a linear regression model as:

$$d\_{ic} = \mathbf{X}\_{ic}\boldsymbol{\mathcal{B}}\_c + e\_{ic} \; E\left(l\_{ic}\right) = \mathbf{X}\_{ic}\boldsymbol{\mathcal{B}}\_c \; e\_{ic} \sim N\left(0, \; 1\right),$$

where Xic corresponds to row i of the design matrix Xc for tumor class c. The link function of the expectation of the liability *Xicβc*with the binary response yic was constructed via a probit model (West, 2003) yielding the following equations:

$$p\_{ic}(y\_{ic} = 1) = \mathsf{O}(X\_{ic}\mathcal{J}\_c) \text{ and } p\_{ic}(y\_{ic} = 1) = 1 - \mathsf{O}(X\_{ic}\mathcal{J}\_c)$$

where *Φ* is the standard normal distribution function. Subject i was classified as having tu‐ mor class c if *pic*(*y*ic=1) was the maximum of the vector pi, containing all *pic*(*y*ic=1) c=1,…, *nc*, where *nc* is the number of tumor classes in the data set.

1 1 )/ 1 )/ 2

where *Pic*contains principle component values for sample i for tumor class c; *βc*is a vector of coefficients estimated using TDS; *nc* is the number of samples in VDS having tumor class *c*;

Following the update of pheromone levels according to equation (2), the PDF is updated ac‐ cording to equation (1) and the process is repeated until some convergence criteria are met. As the PDF is updated, the selected features that perform better will be sampled at higher likelihoods by subsequent artificial ants which, in turn, deposit more "pheromone", thus leading to a positive feedback system similar to the method of communication observed in real ant colonies. Upon convergence the optimal subset of features is select based on the lev‐

*A.4 GCM data set*: The data set contained 198 samples collected from 14 tumor types: BR (breast adenocarcinoma), Pr (prostate adenocarcinoma), LU (lung adenocarcinoma), CO (colorectal adenocarcinoma), LY (lymphoma), BL (bladder transitional cell carcinoma), ML (melanoma), UT (uterine adenocarcinoma), LU (leukemia), RE (renal cell carcinoma), PA (pancreatic adenocarcinoma), OV (ovarian adenocarcinoma), ME (pleural mesothelioma), and CNS (central nervous system). The unedited data set contained the intensity values of 16063 probes generate using Affymetrix high density oligonucleotide microarrays, and cal‐ culated using Affymetrix GENECHIP software (Ramaswamy et al, 2001). Following the thresholding of intensity values to a minimum value of 20 and a maximum value of 16000, a log base 2 transformation was applied to the data set. Genes with the highest expression val‐ ues being less than two times the smallest were removed, leaving 14525 probes for analysis.

*A.5 Results and discussions:* The GCM data set has been a benchmark to compare the perform‐ ance of classification and feature selection algorithms. Table 1 shows the best prediction ac‐ curacies obtained by methods used in this study and several previous studies (GASS (Lin et al., 2006), GA/MLHD (Ooi and Tan, 2003), MAMA (Antonov et al., 2004), and GA/SVM (Liu et al., 2005)) using independent test, performed on the same training and validation data sets originally formed by Ramaswamy et al., 2001 (GCM split), and leave one out cross vali‐ dation (LOOCV). The proposed ACA/LVM yielded substantial increases in accuracies over all other methods, with a 6.5% increase in accuracy over the next best results obtained using the GCM split (Antonov et al., 2004). Furthermore, the ACA/LVM achieved increases of 13.9%, 40%, and 16.6% in accuracy over the FC/LVM, T/LVM, and PT/LVM methods of fea‐

*nc nr acc* = = + - <sup>=</sup> å å **Φ(P βic c Φ(P βic c** (3)

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

181

*nc nr i i*

*c*

and *nr* is the remaining number of samples in VDS.

el of pheromone trail deposited on each feature.

*Δτmc*(*t*)=*accc*

(1−*accc*)

ture selection, respectively.

**6.** The change in pheromone for each tumor class is calculated as:

where accc is the accuracy for tumor type c as calculated using equation (3).

*A.2 Gene Selection*: Filter and wrapper based methods were used to select features to form classifiers for each tumor class. Filter methods selected genes based on ranks determined by the sorted absolute values of fold changes (FC), t-statistics (T), and penalized t-statistics (PT) calculated for each gene for each tumor class. The wrapper method coupled the ACA with LVM (ACA/LVM) such that groups of genes were selected using the ACA and evaluated for performance using LVM.

*A.3 Ant colony optimization*: The general ACA presented in the previous section was used. The prior information,*ηmc* , was assumed as:

$$\eta\_{mc} = \frac{f\_{mc} - \min(f\_c)}{\max(f\_c) - \min(f\_c)} + \frac{t\_{mc} - \min(t\_c)}{\max(t\_c) - \min(t\_c)} + \frac{pt\_{mc} - \min(pt\_c)}{\max(pt\_c) - \min(pt\_c)}}{3}$$

where *f <sup>c</sup>*is a vector of all fold change values for tumor class c; *tc*is a vector of all tstatistic values for tumor class c; and *ptc* is a vector of all penalized t-statistic values for tumor class c. After several trail runs the parameters *α* and *β* were set to 1 and.3 respectively.

The ACA was initialized with all features having an equal baseline level of pheromone used to compute *Pm*(0) for all features. Using the PDF as defined in equation (1), each of *j* artificial ants will select a subset *Sk* of *n* features from the sample space *S* containing all features. The pheromone level of each feature *m* in *Sk* is then updated according to the performance of *Sk* following equation (2).

The procedure can be summarized in the following steps:


$$acc\_c = \frac{\sum\_{i=1}^{nc} \Phi(\mathbf{P\_{ic}} \boldsymbol{\mathfrak{f}\_c}) / nc + \sum\_{i=1}^{nc} 1 - \Phi(\mathbf{P\_{ic}} \boldsymbol{\mathfrak{f}\_c}) / nr}{2} \tag{3}$$

where *Pic*contains principle component values for sample i for tumor class c; *βc*is a vector of coefficients estimated using TDS; *nc* is the number of samples in VDS having tumor class *c*; and *nr* is the remaining number of samples in VDS.

**6.** The change in pheromone for each tumor class is calculated as:

*Δτmc*(*t*)=*accc* (1−*accc*)

where *Φ* is the standard normal distribution function. Subject i was classified as having tu‐ mor class c if *pic*(*y*ic=1) was the maximum of the vector pi, containing all *pic*(*y*ic=1) c=1,…, *nc*,

*A.2 Gene Selection*: Filter and wrapper based methods were used to select features to form classifiers for each tumor class. Filter methods selected genes based on ranks determined by the sorted absolute values of fold changes (FC), t-statistics (T), and penalized t-statistics (PT) calculated for each gene for each tumor class. The wrapper method coupled the ACA with LVM (ACA/LVM) such that groups of genes were selected using the ACA and evaluated for

*A.3 Ant colony optimization*: The general ACA presented in the previous section was used.

where *f <sup>c</sup>*is a vector of all fold change values for tumor class c; *tc*is a vector of all tstatistic values for tumor class c; and *ptc* is a vector of all penalized t-statistic values for tumor class c. After several trail runs the parameters *α* and *β* were set to 1 and.3

The ACA was initialized with all features having an equal baseline level of pheromone used to compute *Pm*(0) for all features. Using the PDF as defined in equation (1), each of *j* artificial ants will select a subset *Sk* of *n* features from the sample space *S* containing all features. The pheromone level of each feature *m* in *Sk* is then updated according to the performance of *Sk*

**2.** Training data is randomly split into two subsets for training (TDS) and validation (VDS) containing ¾ and ¼ of the data, respectively (none of the original validation data

**3.** Using the spectral decomposition of TDS, principle components are computed to allevi‐ ate effects of collinearity and selected for TDS and VDS by removing components with

**4.** Using TDS, a latent variable model is trained for each tumor class, and *p*ic(yic=1) is pre‐

*ptmc* −min(*ptc*) max(*ptc*)−min(*ptc*)

*tmc* −min(*tc*) max(*tc*)−min(*tc*) <sup>+</sup>

3

The procedure can be summarized in the following steps:

**1.** Each ant selects a predetermined number of genes.

(VD) is used at any point in the ACA).

corresponding eigenvalues close to zero.

dicted for every tumor class c for each sample i in VDS.

**5.** The accuracy for each tumor class c is calculated as:

where *nc* is the number of tumor classes in the data set.

180 Ant Colony Optimization - Techniques and Applications

performance using LVM.

*f mc* −min( *f <sup>c</sup>*) max( *<sup>f</sup> <sup>c</sup>*)−min( *<sup>f</sup> <sup>c</sup>*) <sup>+</sup>

*ηmc* =

respectively.

following equation (2).

The prior information,*ηmc* , was assumed as:

where accc is the accuracy for tumor type c as calculated using equation (3).

Following the update of pheromone levels according to equation (2), the PDF is updated ac‐ cording to equation (1) and the process is repeated until some convergence criteria are met. As the PDF is updated, the selected features that perform better will be sampled at higher likelihoods by subsequent artificial ants which, in turn, deposit more "pheromone", thus leading to a positive feedback system similar to the method of communication observed in real ant colonies. Upon convergence the optimal subset of features is select based on the lev‐ el of pheromone trail deposited on each feature.

*A.4 GCM data set*: The data set contained 198 samples collected from 14 tumor types: BR (breast adenocarcinoma), Pr (prostate adenocarcinoma), LU (lung adenocarcinoma), CO (colorectal adenocarcinoma), LY (lymphoma), BL (bladder transitional cell carcinoma), ML (melanoma), UT (uterine adenocarcinoma), LU (leukemia), RE (renal cell carcinoma), PA (pancreatic adenocarcinoma), OV (ovarian adenocarcinoma), ME (pleural mesothelioma), and CNS (central nervous system). The unedited data set contained the intensity values of 16063 probes generate using Affymetrix high density oligonucleotide microarrays, and cal‐ culated using Affymetrix GENECHIP software (Ramaswamy et al, 2001). Following the thresholding of intensity values to a minimum value of 20 and a maximum value of 16000, a log base 2 transformation was applied to the data set. Genes with the highest expression val‐ ues being less than two times the smallest were removed, leaving 14525 probes for analysis.

*A.5 Results and discussions:* The GCM data set has been a benchmark to compare the perform‐ ance of classification and feature selection algorithms. Table 1 shows the best prediction ac‐ curacies obtained by methods used in this study and several previous studies (GASS (Lin et al., 2006), GA/MLHD (Ooi and Tan, 2003), MAMA (Antonov et al., 2004), and GA/SVM (Liu et al., 2005)) using independent test, performed on the same training and validation data sets originally formed by Ramaswamy et al., 2001 (GCM split), and leave one out cross vali‐ dation (LOOCV). The proposed ACA/LVM yielded substantial increases in accuracies over all other methods, with a 6.5% increase in accuracy over the next best results obtained using the GCM split (Antonov et al., 2004). Furthermore, the ACA/LVM achieved increases of 13.9%, 40%, and 16.6% in accuracy over the FC/LVM, T/LVM, and PT/LVM methods of fea‐ ture selection, respectively.


**True\ Predicted**

(50 genes)

(10 genes)

**True\ Predicted** **BR PR LU CO LY BL ML UT LE RE PA OV ME CNS**

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

183

1 6 4 6 6 7 2 2 6 3 1 2 4 4 43/54

1 6 4 6 6 7 2 2 6 3 1 2 4 4 42/54

**Table 4.** Confusion matrix for best predictions obtained for GCM data set using genes selected by the penalized t-test

**Table 3.** Confusion matrix for best predictions obtained for the GCM data set using genes selected by the fold change

**BR PR LU CO LY BL ML UT LE RE PA OV ME CNS**

BR 0 3 1 4 PR 1 5 6 LU 4 4 CO 4 4 LY 6 6 BL 1 2 3 ML 2 2 UT 2 2 LE 6 6 RE 2 1 3 PA 2 1 0 3 OV 1 2 1 4 ME 3 3 CNS 4 4

BR 0 3 1 4 PR 1 5 6 LU 3 1 4 CO 4 4 LY 6 6 BL 1 2 3 ML 2 2 UT 2 2 LE 6 6 RE 2 1 3 PA 1 1 1 3 OV 1 3 1 4 ME 3 3 CNS 4 4

a Split used by Ramaswamy et al 2001; bLeave one out cross validation; c Number of genes selected prior to the imple‐ mentation of feature selection algorithm; dWeighted average of scaled fold change, t-test, and penalized t-test values.

**Table 1.** Accuracy (%) of tumor class predictions using ant colony algorithm (ACA) and several previously published methods.

Due to its poor performance, the confusion matrix of predictions using T/LVM is not includ‐ ed, but matrices for the predictions obtained by the ACA/LVM, FC/LVM, and PT/LVM us‐ ing the GCM split can be found in Tables 2-4. These tables show that the ACA/LVM performs as good or better than the rank based methods for every tumor type. Additionally the ACA/LVM correctly predicted 50% of the BR samples, a tumor class that has traditional‐ ly yielded very poor results (Bagirov et al., 2003; Ramaswamy et al., 2001). The ACA/LVM also achieved 100% prediction accuracy for 10 of the 14 tumor classes, as compared to only 7 and 8 when using FC/LVM or PT/LVM, respectively.


**Table 2.** Confusion matrix for predictions obtained for the GCM data set using genes selected by the ant colony algorithm.


**GCM data set**

mentation of feature selection algorithm; dWeighted average of scaled fold change, t-test, and penalized t-test values.

**Table 1.** Accuracy (%) of tumor class predictions using ant colony algorithm (ACA) and several previously published

Due to its poor performance, the confusion matrix of predictions using T/LVM is not includ‐ ed, but matrices for the predictions obtained by the ACA/LVM, FC/LVM, and PT/LVM us‐ ing the GCM split can be found in Tables 2-4. These tables show that the ACA/LVM performs as good or better than the rank based methods for every tumor type. Additionally the ACA/LVM correctly predicted 50% of the BR samples, a tumor class that has traditional‐ ly yielded very poor results (Bagirov et al., 2003; Ramaswamy et al., 2001). The ACA/LVM also achieved 100% prediction accuracy for 10 of the 14 tumor classes, as compared to only 7

**BR PR LU CO LY BL ML UT LE RE PA OV ME CNS**

1 6 4 6 6 7 2 2 6 3 1 2 4 4 49/54

**Table 2.** Confusion matrix for predictions obtained for the GCM data set using genes selected by the ant colony algorithm.

BR 2 1 1 4 PR 1 5 6 LU 4 4 CO 4 4 LY 6 6 BL 1 2 3 ML 2 2 UT 2 2 LE 6 6 RE 3 3 PA 1 2 3 OV 4 4 ME 3 3 CNS 4 4

FC/LVM(14525) 79.6 74.8 \_\_\_\_ T/LVM(14525) 64.8 \_\_\_\_ \_\_\_\_ PT/LVM(14525) 77.8 74.4 \_\_\_\_ AVGd/LVM(14525) 79.6 74.8 \_\_\_\_ GASS(1000) 81.5 \_\_\_\_ 81.3 GA/MLHD(1000) 76 \_\_\_\_ 79.8 MAMA 85.2 \_\_\_\_ \_\_\_\_\_ GA/SVM(1000) \_\_\_ \_\_\_\_ 81

Split used by Ramaswamy et al 2001; bLeave one out cross validation; c

and 8 when using FC/LVM or PT/LVM, respectively.

) 90.7 84.8 \_\_\_\_

ACA/LVM(14525c

182 Ant Colony Optimization - Techniques and Applications

a

methods.

**True\ Predicted** GCM splita Replicated splits LOOCVb

Number of genes selected prior to the imple‐

**Table 3.** Confusion matrix for best predictions obtained for the GCM data set using genes selected by the fold change (50 genes)


**Table 4.** Confusion matrix for best predictions obtained for GCM data set using genes selected by the penalized t-test (10 genes)

To further evaluate performance, each of the feature selection algorithms was tested using four additional random splits of the data. The best classification accuracies obtained for each algorithm can be found in Table 5. The ACA/LVM algorithm yielded the best prediction ac‐ curacies for all replicates, with increases in accuracies ranging from 6.7% to 14% over the best accuracies obtained by filter methods. When looking at the three filter methods it can be seen that the best method varied depending on the replication. These findings are in agree‐ ment with Jefferey et al. (2006).

**BR PR LU CO LY BL ML UT LE RE PA OV ME CNS**

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

185

ACA 3.4 4.8 2 7.8 6.6 19.6 4.6 7.6 3.2 16 14.6 17.2 5 5.6 FC 18 18 18 18 18 18 18 18 18 18 18 18 18 18 PT 14 14 14 14 14 14 14 14 14 14 14 14 14 14 Averagea 18 18 18 18 18 18 18 18 18 18 18 18 18 18

Weighted average of scaled fold change (FC), t-test, and penalized t-test (PT) values

making the ACA a versatile feature selection tool.

**Table 6.** Number of genes selected for each tumor type using ACA and other feature selection methods.

The superiority of the ACA/LVM when compared to models using GA indicates the ACA's utility, as compared to other optimization methods, when working with high dimension da‐ ta sets. The ACA's ability to incorporate prior information in the optimization process pro‐ vides several advantages over other optimization algorithms when dealing with large numbers of features. The inclusion of prior information in the pheromone function focuses the selection process on genes that should yield better results without the need for an explic‐ it truncation of the data, which was needed to achieve good results with the GA (Hong and Cho, 2006; Lin et al., 2006; Liu et al., 2005; Ooi and Tan et al., 2003; Peng et al., 2003). Trunca‐ tion of large numbers of genes could a priori eliminate genes from consideration that, though they may not have high predictive ability alone, could contribute to the predictive power of an ensemble of genes. Additionally, depending on the method of truncation, the reduced gene list could be highly redundant (Lin et al., 2006; Shen et al., 2006), further re‐ ducing the informativeness of pre-selected genes. Conversely, when removing a small num‐ ber of features in a large data set, the truncated data set may be too large for efficient convergence of the algorithm (Lin et al., 2006). Additionally, the inclusion of prior informa‐ tion allows the ACA to be coupled with many other types of feature selection methods,

For LU tumors, the ACA identified two genes capable of classifying LU tumor samples with 100%, in each of the five replicates. The selected genes, SP-B and SP-A, both encode pulmo‐ nary surfactant proteins which are necessary for lung function. Another tumor class, with which the ACA was able to select a small number of highly predictive genes, was CNS. As with the LU tumor type, the genes selected by the ACA were very consistent from replica‐ tion to replication. The gene encoding for APCL protein had the highest pheromone levels in all five replicates and was the only gene required to achieve 100% accuracy in replicate five. APCL protein is a homologue of APC, a known tumor suppressor that interacts with microtubules during mitosis (Akiyama and Kawasaki, 2006). The gene encoding MAP1B, a protein found to be important in synaptic function of cortical neurons, was also identified as being highly predictive of CNS tumor types. Several other genes selected by the ACA, found

in *supplemental materials*, were identified in a previous study (Antonov et al., 2004).

In contrast to the LU and CNS tumor types, BR samples were consistently predicted with low accuracies. These findings are in agreement with previous results (Bagirov et al., 2003; Ramaswamy et al., 2001). Unlike the gene list obtained for BR and CNS tumor types, the

a


a Split used by Ramaswamy et al 2001; bWeighted average of scaled fold change (FC),

t-test (PT), and penalized t-test values (PT).

**Table 5.** Classification accuracies using several feature selection methods

Due to a lack of any good criterion for determining an objective cut-off value for the rank based methods, several values were used and evaluated. Since the use of fewer fea‐ tures is desirable from a biological standpoint, an upper limit of 50 genes per tumor class was imposed on all methods. Table 6 shows the number of genes needed for each tumor type to achieve the best results, averaged across all replicates. It can be seen that, for 10 of the 14 tumor classes, the ACA/LVM selects fewer genes than the rank based methods.

The performance of the ACA/LVM model was superior, not only to the filter based methods used in this study, but also several reported results using the GCM data set. The ACA/LVM consistently yielded superior accuracies using fewer genes than the filter based methods, for which ranks varied with each replication. The breaks in pheromone levels observed with the most predictive genes also provided more objective selection criteria for identifying top fea‐ tures, unlike the filter methods in which truncation points were somewhat arbitrary. The ob‐ jective selection criteria and robustness of the ACA, within the confines of the GCM data set, make it a superior method for clinical applications, as it could enable a single procedure to be effectively applied to varied applications. The use of filter based methods in such scenar‐ ios would require different combinations of truncation points and scoring methods for each data set, a highly impractical endeavor.


**Table 6.** Number of genes selected for each tumor type using ACA and other feature selection methods.

To further evaluate performance, each of the feature selection algorithms was tested using four additional random splits of the data. The best classification accuracies obtained for each algorithm can be found in Table 5. The ACA/LVM algorithm yielded the best prediction ac‐ curacies for all replicates, with increases in accuracies ranging from 6.7% to 14% over the best accuracies obtained by filter methods. When looking at the three filter methods it can be seen that the best method varied depending on the replication. These findings are in agree‐

**Replication 1 2 3 4 5**

ACA/LVM 90.7 83.3 79.6 81.5 88.9

FC/LVM 79.6 77.8 68.5 72.2 75.9

PT/LVM 77.8 77.8 66.7 68.5 81.5

AVGb/LVM 79.6 70.4 70.4 70.4 83.3

Due to a lack of any good criterion for determining an objective cut-off value for the rank based methods, several values were used and evaluated. Since the use of fewer fea‐ tures is desirable from a biological standpoint, an upper limit of 50 genes per tumor class was imposed on all methods. Table 6 shows the number of genes needed for each tumor type to achieve the best results, averaged across all replicates. It can be seen that, for 10 of the 14 tumor classes, the ACA/LVM selects fewer genes than the rank based

The performance of the ACA/LVM model was superior, not only to the filter based methods used in this study, but also several reported results using the GCM data set. The ACA/LVM consistently yielded superior accuracies using fewer genes than the filter based methods, for which ranks varied with each replication. The breaks in pheromone levels observed with the most predictive genes also provided more objective selection criteria for identifying top fea‐ tures, unlike the filter methods in which truncation points were somewhat arbitrary. The ob‐ jective selection criteria and robustness of the ACA, within the confines of the GCM data set, make it a superior method for clinical applications, as it could enable a single procedure to be effectively applied to varied applications. The use of filter based methods in such scenar‐ ios would require different combinations of truncation points and scoring methods for each

Split used by Ramaswamy et al 2001; bWeighted average of scaled fold change (FC),

**Table 5.** Classification accuracies using several feature selection methods

ment with Jefferey et al. (2006).

184 Ant Colony Optimization - Techniques and Applications

t-test (PT), and penalized t-test values (PT).

data set, a highly impractical endeavor.

a

methods.

The superiority of the ACA/LVM when compared to models using GA indicates the ACA's utility, as compared to other optimization methods, when working with high dimension da‐ ta sets. The ACA's ability to incorporate prior information in the optimization process pro‐ vides several advantages over other optimization algorithms when dealing with large numbers of features. The inclusion of prior information in the pheromone function focuses the selection process on genes that should yield better results without the need for an explic‐ it truncation of the data, which was needed to achieve good results with the GA (Hong and Cho, 2006; Lin et al., 2006; Liu et al., 2005; Ooi and Tan et al., 2003; Peng et al., 2003). Trunca‐ tion of large numbers of genes could a priori eliminate genes from consideration that, though they may not have high predictive ability alone, could contribute to the predictive power of an ensemble of genes. Additionally, depending on the method of truncation, the reduced gene list could be highly redundant (Lin et al., 2006; Shen et al., 2006), further re‐ ducing the informativeness of pre-selected genes. Conversely, when removing a small num‐ ber of features in a large data set, the truncated data set may be too large for efficient convergence of the algorithm (Lin et al., 2006). Additionally, the inclusion of prior informa‐ tion allows the ACA to be coupled with many other types of feature selection methods, making the ACA a versatile feature selection tool.

For LU tumors, the ACA identified two genes capable of classifying LU tumor samples with 100%, in each of the five replicates. The selected genes, SP-B and SP-A, both encode pulmo‐ nary surfactant proteins which are necessary for lung function. Another tumor class, with which the ACA was able to select a small number of highly predictive genes, was CNS. As with the LU tumor type, the genes selected by the ACA were very consistent from replica‐ tion to replication. The gene encoding for APCL protein had the highest pheromone levels in all five replicates and was the only gene required to achieve 100% accuracy in replicate five. APCL protein is a homologue of APC, a known tumor suppressor that interacts with microtubules during mitosis (Akiyama and Kawasaki, 2006). The gene encoding MAP1B, a protein found to be important in synaptic function of cortical neurons, was also identified as being highly predictive of CNS tumor types. Several other genes selected by the ACA, found in *supplemental materials*, were identified in a previous study (Antonov et al., 2004).

In contrast to the LU and CNS tumor types, BR samples were consistently predicted with low accuracies. These findings are in agreement with previous results (Bagirov et al., 2003; Ramaswamy et al., 2001). Unlike the gene list obtained for BR and CNS tumor types, the gene lists for BR tumors were highly variable, suggesting potentially high heterogeneity in these tumor samples. Despite dissimilarities between the genes selected across replications, the ACA did identify SEPT9 as being highly predictive in four of the five replicates. The pro‐ tein encoded by this gene has been shown to be involved in mitosis of mammary epithelial cells (Nagata et al., 2003) and has been associated with both ovarian and breast neoplasia (Scott et al., 2006). The identification of this gene by the ACA demonstrates its ability to identify biologically relevant features in challenging data sets.

*yi* ={

where *P*<sup>i</sup>

groups.

equations:

yi={

1 if *lori* ≥0 0 if *lori* <0

The log odds ratio *lor*<sup>i</sup>

= probability (yi

The link function of the log odds ratio *X<sup>i</sup>*

(*y*i=0)= <sup>1</sup>

1 + exp(*X<sup>i</sup>*

*pi*

yielding the following relationships:

*β*)

*β*)

*<sup>β</sup>*) <sup>≥</sup>0.5

*<sup>β</sup>*) <0.5

<sup>1</sup> if exp(*X*<sup>i</sup>

<sup>0</sup> if exp(*X*<sup>i</sup>

1 + exp(*X*<sup>i</sup>

1 + exp(*X*<sup>i</sup>

is modeled as:

*lori* =ln( *pi*

1− *pi*

*<sup>β</sup>*) and *pi*

*B.2 Marginal effects model*: The genotype and haplotype association methods were imple‐ mented using R functions developed by [Gonzalez et al., 2007; Sinnwell and Schaid, 2005]. The haplotype analysis was implemented using a sliding window approach which utilizes a window of *k* SNP in width sliding across the genome *h* SNP at a time. Individual SNP scores

*B.3 Ant colony algorithm*: While the algorithm, in the aforementioned form can be used to subjectively identify markers, it is not well suited for the calculation of permutation p-val‐ ues. When updating the pheromone function, as previously described in equation (2), the fi‐ nal pheromone levels are relative not only to prediction accuracy, but the number of times a SNP marker is selected. As a result, the amount of pheromone deposited on a feature de‐ pends greatly on the amount of pheromone deposited on all other SNP markers and can vary wildly from permutation to permutation. One obvious solution to this problem is to use the average accuracy of all *Sk* containing genotypes for SNP *m*; however, this approach substantially reduces the ACA's ability to efficiently burn in on good solutions, an attribute

were determined as the maximum average of all haplotypes containing a given SNP.

needed to detect unknown gene interactions in high-dimension data sets.

To overcome these limitations, a two-layer pheromone function was developed:

)= *X<sup>i</sup>*

types formed from the selected SNP. Groups of SNP markers with less than two correspond‐ ing observations were discarded, and analysis was conducted on all remaining marker

*β* + *ei* (4)

Ant Colony Algorithm with Applications in the Field of Genomics

*β*with the binary response yi gives the following

*<sup>β</sup>*) (5)

http://dx.doi.org/10.5772/52051

187

*β*)

1 + exp(*X<sup>i</sup>*

= 1) and X is a matrix containing indicator variables for the haplo‐

(*y*i=1)= exp(*X<sup>i</sup>*

## **2.2. The use of the ant colony algorithm for the detection of marker associations in the presence of gene interactions**

With the advent of high-throughput, cost effective genotyping platforms, there has been much focus on the use of high-density single nucleotide polymorphism (SNP) genotyp‐ ing to identify causative mutations for traits of interest, and while putative mutations have been identified for several traits, these studies tend to focus on SNP with large marginal effects [Hugot et al., 2001; Woon et al., 2007]. However, several studies have found that gene interactions may play important roles in many complex traits [Coutinho et al., 2007; Barendse et al., 2007]. Given the high density of SNP maker maps, examin‐ ing all possible interactions is seldom possible computationally. As a result, studies ex‐ amining gene interactions tend to focus on a small number of SNP, previously identified as having strong marginal associations. Using an exhaustive search of all two-way inter‐ actions, Marchini et al. achieved greater power to detect causative mutations than when estimating only marginal effects. Due to the high computational cost of this approach, a two-stage model was proposed, in which SNP were selected in the first stage based on marginal effects and then tested for interactions in the subsequent stage [Marchini et al., 2005]. This approach could, however, result in the failure to detect important regions of the genome in the first stage of the model. As such, there is a need for methodologies capable of identifying important genomic regions in the presence of potential gene inter‐ actions when large numbers of markers are genotyped.

One approach would be to view the identification of groups of interacting SNP as an optimi‐ zation problem, for which several algorithms have been developed. These algorithms are designed to search large sample spaces for globally optimal solutions and have been applied to a wide range of problems [Shymygelska and Hoos, 2005; Ding et al., 2005]. Through the evaluation of groups of loci efficiently selected from different regions of the genome, optimi‐ zation algorithms should be able to account for potential interactions.

In this section, a modified ACA, enabling the use of permutation testing for global signifi‐ cance, was combined with logistic regression and implemented on a simulated binary trait under the influence of interacting genes. The performance of the ACA was evaluated and compared to models accounting for only marginal effects.

*B.1 Logistic regression*: Groups of SNP markers were evaluated based in haplotype genotype effects estimated as log odds ratios (*lor*) using logistic regression (LR). The relationship be‐ tween the *lor* and the binary response can be expressed as:

$$y\_i = \begin{cases} 1 & \text{if } \operatorname{lor}\_i \ge 0 \\ 0 & \text{if } \operatorname{lor}\_i < 0 \end{cases}$$

gene lists for BR tumors were highly variable, suggesting potentially high heterogeneity in these tumor samples. Despite dissimilarities between the genes selected across replications, the ACA did identify SEPT9 as being highly predictive in four of the five replicates. The pro‐ tein encoded by this gene has been shown to be involved in mitosis of mammary epithelial cells (Nagata et al., 2003) and has been associated with both ovarian and breast neoplasia (Scott et al., 2006). The identification of this gene by the ACA demonstrates its ability to

**2.2. The use of the ant colony algorithm for the detection of marker associations in the**

With the advent of high-throughput, cost effective genotyping platforms, there has been much focus on the use of high-density single nucleotide polymorphism (SNP) genotyp‐ ing to identify causative mutations for traits of interest, and while putative mutations have been identified for several traits, these studies tend to focus on SNP with large marginal effects [Hugot et al., 2001; Woon et al., 2007]. However, several studies have found that gene interactions may play important roles in many complex traits [Coutinho et al., 2007; Barendse et al., 2007]. Given the high density of SNP maker maps, examin‐ ing all possible interactions is seldom possible computationally. As a result, studies ex‐ amining gene interactions tend to focus on a small number of SNP, previously identified as having strong marginal associations. Using an exhaustive search of all two-way inter‐ actions, Marchini et al. achieved greater power to detect causative mutations than when estimating only marginal effects. Due to the high computational cost of this approach, a two-stage model was proposed, in which SNP were selected in the first stage based on marginal effects and then tested for interactions in the subsequent stage [Marchini et al., 2005]. This approach could, however, result in the failure to detect important regions of the genome in the first stage of the model. As such, there is a need for methodologies capable of identifying important genomic regions in the presence of potential gene inter‐

One approach would be to view the identification of groups of interacting SNP as an optimi‐ zation problem, for which several algorithms have been developed. These algorithms are designed to search large sample spaces for globally optimal solutions and have been applied to a wide range of problems [Shymygelska and Hoos, 2005; Ding et al., 2005]. Through the evaluation of groups of loci efficiently selected from different regions of the genome, optimi‐

In this section, a modified ACA, enabling the use of permutation testing for global signifi‐ cance, was combined with logistic regression and implemented on a simulated binary trait under the influence of interacting genes. The performance of the ACA was evaluated and

*B.1 Logistic regression*: Groups of SNP markers were evaluated based in haplotype genotype effects estimated as log odds ratios (*lor*) using logistic regression (LR). The relationship be‐

identify biologically relevant features in challenging data sets.

actions when large numbers of markers are genotyped.

compared to models accounting for only marginal effects.

tween the *lor* and the binary response can be expressed as:

zation algorithms should be able to account for potential interactions.

**presence of gene interactions**

186 Ant Colony Optimization - Techniques and Applications

The log odds ratio *lor*<sup>i</sup> is modeled as:

$$\text{Mor}\_i = \ln(\frac{p\_i}{1 - p\_i}) = \mathbf{X}\_i \boldsymbol{\mathcal{B}} + \boldsymbol{e}\_i \tag{4}$$

where *P*<sup>i</sup> = probability (yi = 1) and X is a matrix containing indicator variables for the haplo‐ types formed from the selected SNP. Groups of SNP markers with less than two correspond‐ ing observations were discarded, and analysis was conducted on all remaining marker groups.

The link function of the log odds ratio *X<sup>i</sup> β*with the binary response yi gives the following equations:

$$p\_i(y\_i = 0) = \frac{1}{1 + \exp(X\_i \beta)} \text{ and } p\_i(y\_i = 1) = \frac{\exp(X\_i \beta)}{1 + \exp(X\_i \beta)}\tag{5}$$

yielding the following relationships:

$$\mathbf{y}\_{i} = \begin{cases} 1 & \text{if } \frac{\exp(\mathbf{X}\_{i}\boldsymbol{\beta})}{1 + \exp(\mathbf{X}\_{i}\boldsymbol{\beta})} \succeq 0.5 \\\\ 0 & \text{if } \frac{\exp(\mathbf{X}\_{i}\boldsymbol{\beta})}{1 + \exp(\mathbf{X}\_{i}\boldsymbol{\beta})} < 0.5 \end{cases}$$

*B.2 Marginal effects model*: The genotype and haplotype association methods were imple‐ mented using R functions developed by [Gonzalez et al., 2007; Sinnwell and Schaid, 2005]. The haplotype analysis was implemented using a sliding window approach which utilizes a window of *k* SNP in width sliding across the genome *h* SNP at a time. Individual SNP scores were determined as the maximum average of all haplotypes containing a given SNP.

*B.3 Ant colony algorithm*: While the algorithm, in the aforementioned form can be used to subjectively identify markers, it is not well suited for the calculation of permutation p-val‐ ues. When updating the pheromone function, as previously described in equation (2), the fi‐ nal pheromone levels are relative not only to prediction accuracy, but the number of times a SNP marker is selected. As a result, the amount of pheromone deposited on a feature de‐ pends greatly on the amount of pheromone deposited on all other SNP markers and can vary wildly from permutation to permutation. One obvious solution to this problem is to use the average accuracy of all *Sk* containing genotypes for SNP *m*; however, this approach substantially reduces the ACA's ability to efficiently burn in on good solutions, an attribute needed to detect unknown gene interactions in high-dimension data sets.

To overcome these limitations, a two-layer pheromone function was developed:

$$P\_m(t) = \frac{\tau\_m(t)^\alpha \tau \mathcal{D}\_m(t)^{\alpha 2} \eta\_m \beta}{\sum\_{m=1}^{n^f} \tau\_m(t)^\alpha \tau \mathcal{D}\_m(t)^{\alpha 2} \eta\_m \beta} \tag{6}$$

**Scenario 1 Scenario 2** AB aB Ab ab AB aB Ab ab

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

189

AB 1 1 1 1 1 1 1 1 aB 1 1 1 1 1 1 1 1 Ab 1 1 1 1 1 1 1 1 Ab 1 1 1 15 1 1 1 10

The loci of the causative mutations were selected at random; with the frequencies of the causa‐ tive mutations being.58 and.6. Although these frequencies might be considered high, it was necessary to restrict selection to SNP with mutant allele frequencies greater than.5. This was done to insure a reasonable simulated disease incidence of 15%. A plot illustrating the LD of all SNP with the two causative mutations is shown in Fig (1). The plot shows a large peak of high LD with rs2049736 (SNP 409), while the peak of high LD with rs28953468 (SNP 2041) is substan‐

**Figure 1.** Plots of each marker's linkage disequilibrium (LD) with the two causative mutations. The light grey line rep‐ resents LD with the causative mutation located at position 409. The black line represents LD with the causative muta‐

tion located at position 2041.

tially narrower, and is preceded by a plateau of SNP in moderate LD with rs28953468.

**Table 7.** Relative risk for simulated trait (relative to the aa/bb genotype)

where *τm*(*t*) is the first pheromone layer updated using the sum of accuracies for all *Sk* con‐ taining SNP *m*; *τ*2*m*(*t*)is the second pheromone layer updated using the average accuracy of all *Sk* containing genotypes for SNP *m*; and*ηm*, *α*, *β*are as previously described. For the cur‐ rent study, *α*and *α*2 were set to 1, *β*was set to.3 and the prior information (*ηm*) was the pre‐ diction the accuracy of SNP marker *m*, obtained using logistic regression on genotypes.

The pheromone for *τm*(*t*) was updated using equation (2) and *τ*2*m*(*t*) was updated using the following equation:

$$
\tau \mathcal{Z}\_m(t+1) = \left[t \, \, ^\ast \tau\_m \mathcal{Z}(t) + \Delta \tau\_m \mathcal{Z}(t)\right] / \left(t + ns\right) \tag{7}
$$

where *t* is the iteration number; *Δτm*2(*t*)is the change in pheromone level for feature *m* based on the sum of accuracy of all *Sk* containing genotypes for SNP *m*, and is set to zero if feature *m* was not selected by any of the artificial ants; and ns is the number of times SNP *m* was selected at iteration *t*. Permutation p-values were calculated using *τ*2*m*(*t*) only.

The procedure can be summarized in the following steps:


$$
\omega\_{\text{pero}\,\text{moone}} = \text{acc}^{\text{(1-acc)}} \tag{8}
$$


*B.4 Data simulation*: Genotype data on 90 unrelated individuals from the Japanese and Han Chinese populations were downloaded from the HapMap ECODE project website. Each simulation scenario was replicated five times using two 500 Kbp regions on chro‐ mosome 2, comprising 2047 polymorphic SNP. All SNP haplotypes were assumed to be known without error. The binary disease trait was simulated under a two locus epistatic model as seen in Table 7.


**Table 7.** Relative risk for simulated trait (relative to the aa/bb genotype)

2

 h

 ab

() 2 () *m mm*

*t t*

*m m mm*

=

where *τm*(*t*) is the first pheromone layer updated using the sum of accuracies for all *Sk* con‐ taining SNP *m*; *τ*2*m*(*t*)is the second pheromone layer updated using the average accuracy of all *Sk* containing genotypes for SNP *m*; and*ηm*, *α*, *β*are as previously described. For the cur‐ rent study, *α*and *α*2 were set to 1, *β*was set to.3 and the prior information (*ηm*) was the pre‐ diction the accuracy of SNP marker *m*, obtained using logistic regression on genotypes.

The pheromone for *τm*(*t*) was updated using equation (2) and *τ*2*m*(*t*) was updated using the

where *t* is the iteration number; *Δτm*2(*t*)is the change in pheromone level for feature *m* based on the sum of accuracy of all *Sk* containing genotypes for SNP *m*, and is set to zero if feature *m* was not selected by any of the artificial ants; and ns is the number of times SNP *m* was

**2.** Using the selected SNP markers, accuracies are computed using logistic regression on

**1.** The change in pheromone at time *t* is then calculated using equations (2) and (7).

**2.** Following the update of pheromone levels according to equations (2) and (7), the PDF is updated according to equation (6) and the process is repeated until pheromone levels

*B.4 Data simulation*: Genotype data on 90 unrelated individuals from the Japanese and Han Chinese populations were downloaded from the HapMap ECODE project website. Each simulation scenario was replicated five times using two 500 Kbp regions on chro‐ mosome 2, comprising 2047 polymorphic SNP. All SNP haplotypes were assumed to be known without error. The binary disease trait was simulated under a two locus epistatic

(1 ) *acc*

*<sup>k</sup> pheromone acc* - <sup>=</sup> (8)

2 ( 1) [ \* 2( ) 2( )] / ( ) *m mm*

tt

selected at iteration *t*. Permutation p-values were calculated using *τ*2*m*(*t*) only.

a

tt

1 () 2 () ( )

a

tt

*t t P t*

*m nf*

following equation:

t

188 Ant Colony Optimization - Techniques and Applications

haplotypes or genotypes.

have converged.

model as seen in Table 7.

The procedure can be summarized in the following steps:

**1.** Each ant selects a predetermined number of SNP markers.

**3.** The pheromone for each selected group of SNP, *Sk* , is calculated as:

=

2

 h

å (6)

*t t t t t ns* + = +D + (7)

 ab

> The loci of the causative mutations were selected at random; with the frequencies of the causa‐ tive mutations being.58 and.6. Although these frequencies might be considered high, it was necessary to restrict selection to SNP with mutant allele frequencies greater than.5. This was done to insure a reasonable simulated disease incidence of 15%. A plot illustrating the LD of all SNP with the two causative mutations is shown in Fig (1). The plot shows a large peak of high LD with rs2049736 (SNP 409), while the peak of high LD with rs28953468 (SNP 2041) is substan‐ tially narrower, and is preceded by a plateau of SNP in moderate LD with rs28953468.

**Figure 1.** Plots of each marker's linkage disequilibrium (LD) with the two causative mutations. The light grey line rep‐ resents LD with the causative mutation located at position 409. The black line represents LD with the causative muta‐ tion located at position 2041.

Permutation testing was used to access global significance for all models used in the study. Statuses were randomly shuffled amongst subjects, with haplotype effects, genotype effects and association p-values re-estimated for each new configuration of the response variables. The largest estimated haplotype/genotype effect or the smallest haplotype/genotype associa‐ tion p-value from each permutation was saved to form an empirical distribution used for calculation of p-values. One hundred permutations were performed, yielding p-values accu‐ rate to 1%. Power was calculated as the proportion of times a given method identified at least one SNP marker in high LD (r2 ≥.80) with a causative mutation.

*B.5 Results and discussions*: Estimates of power for the three methods can be found in Table 8. Methods employing the ACA showed substantial increases in power when compared to the methods accounting for only marginal effects. Due to the fact that the trait was simulated under a dominance model, analysis of genotypes yielded superior results when compared to haplotype analysis. Despite the inherent advantage of genotype analysis using a dominance model, the ACA using haplotypes (ACA/H) still showed greater power than RG/D in both scenarios. For scenario 2, all models showed a reduction in power; however, the superiority of the ACA methodologies remained constant, with the ACA using LG on genotypes assum‐ ing a dominance model (ACA/G/D) yielding 66.7% increase in power for both scenarios when compared to the next best method, RG/D.

(a) (b)

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

191

(c)

**Figure 2.** Association plots of SNP markers for the simulated trait under scenario 1. Plots were obtained using 2 SNP haplotypes analyzed by a. SW/LR and b. ACA/LR. Vertical lines represent the position of the two causative mutations,

(a) (b)

(c)

**Figure 3.** Association plots of SNP markers for the simulated trait under scenario 2. Plots were obtained using 3 SNP haplotypes analyzed by a. SW/LR, b. ACA/LR, and c. RG. Vertical lines represent the position of the two causative mu‐

tations, and horizontal lines represent the threshold at which associations are significant at α=.05.

and horizontal lines represent the threshold at which associations are significant at α=. 05


<sup>a</sup> Power was calculated as the proportion of times at least one SNP in high linkage disequilibrium (>.8) with a causative mutations was detected by the model at α=.05 for genome-wide significance

#### **Table 8.** Power calculationsa .

Plots of the associative effects, obtained using SW/H, ACA/G/D, and RG/D, are shown in Fig. (2) and (3). When compared to the LD plot (Fig. (1)) all methods show good correspond‐ ence for scenario 1, though only the ACA/G/D was able to identify markers for both causa‐ tive mutations in all replicates. In scenario 2, where the genetic effect was greatly reduced, plots of associative effects tended to be noisier for all models, with the ACA/G/D again showing superior performance, identifying several SNP markers having only moderate LD with causative mutation rs28953468.

Permutation testing was used to access global significance for all models used in the study. Statuses were randomly shuffled amongst subjects, with haplotype effects, genotype effects and association p-values re-estimated for each new configuration of the response variables. The largest estimated haplotype/genotype effect or the smallest haplotype/genotype associa‐ tion p-value from each permutation was saved to form an empirical distribution used for calculation of p-values. One hundred permutations were performed, yielding p-values accu‐ rate to 1%. Power was calculated as the proportion of times a given method identified at

*B.5 Results and discussions*: Estimates of power for the three methods can be found in Table 8. Methods employing the ACA showed substantial increases in power when compared to the methods accounting for only marginal effects. Due to the fact that the trait was simulated under a dominance model, analysis of genotypes yielded superior results when compared to haplotype analysis. Despite the inherent advantage of genotype analysis using a dominance model, the ACA using haplotypes (ACA/H) still showed greater power than RG/D in both scenarios. For scenario 2, all models showed a reduction in power; however, the superiority of the ACA methodologies remained constant, with the ACA using LG on genotypes assum‐ ing a dominance model (ACA/G/D) yielding 66.7% increase in power for both scenarios

ACA/G/D \_\_\_ 1.00 0.90 \_\_\_ 0.50 0.40 ACA/G/C \_\_\_ 0.70 0.80 \_\_\_ 0.40 0.40 ACA/HAP \_\_\_ 0.60 0.70 \_\_\_ 0.50 0.40 RG/D 0.60 \_\_\_ \_\_\_ 0.30 \_\_\_ \_\_\_ RG/C 0.30 \_\_\_ \_\_\_ 0.30 \_\_\_ \_\_\_ SW/HAP \_\_\_ 0.10 0.20 \_\_\_ 0.00 0.00

<sup>a</sup> Power was calculated as the proportion of times at least one SNP in high linkage disequilibrium (>.8) with a causative

Plots of the associative effects, obtained using SW/H, ACA/G/D, and RG/D, are shown in Fig. (2) and (3). When compared to the LD plot (Fig. (1)) all methods show good correspond‐ ence for scenario 1, though only the ACA/G/D was able to identify markers for both causa‐ tive mutations in all replicates. In scenario 2, where the genetic effect was greatly reduced, plots of associative effects tended to be noisier for all models, with the ACA/G/D again showing superior performance, identifying several SNP markers having only moderate LD

≥.80) with a causative mutation.

**Scenario 1 Scenario 2** 1 locus 2 locus 3 locus 1 locus 2 locus 3 locus

least one SNP marker in high LD (r2

190 Ant Colony Optimization - Techniques and Applications

when compared to the next best method, RG/D.

mutations was detected by the model at α=.05 for genome-wide significance

.

with causative mutation rs28953468.

**Table 8.** Power calculationsa

**Figure 2.** Association plots of SNP markers for the simulated trait under scenario 1. Plots were obtained using 2 SNP haplotypes analyzed by a. SW/LR and b. ACA/LR. Vertical lines represent the position of the two causative mutations, and horizontal lines represent the threshold at which associations are significant at α=. 05

**Figure 3.** Association plots of SNP markers for the simulated trait under scenario 2. Plots were obtained using 3 SNP haplotypes analyzed by a. SW/LR, b. ACA/LR, and c. RG. Vertical lines represent the position of the two causative mu‐ tations, and horizontal lines represent the threshold at which associations are significant at α=.05.

To determine the effectiveness of the permutation on pheromone levels, the cumulative dis‐ tribution, based on LD with causative mutations, of SNP identified as being significantly as‐ sociated with simulated trait by ACA/G/D and RG/D were plotted and can be found in Fig. (4). Despite similarities in the average number of SNP identified by ACA/G/D (15.4) and RG/D (22), the distributions of these SNP, differed substantially. In contrast to RG/D, the ACA/G/D identified a large number of SNP having LD between.35-.45. These SNP corre‐ sponded to the broad plateau of SNP in LD with SNP 2041. Unlike RG/D, the ACA/G/D also identified several SNP (5.19%) having less than.10 LD with either of the causative mutations, an unexpected result given the strict family-wise significance thresholds (α=0.05) imposed on all models. Surprisingly, both methodologies identified a large number of SNP having LD of approximately ~.2. Upon closer examination it was found that these SNP had LD of ~. 2 with both causative mutations, likely artifacts of the data resulting from the relatively small sample size. The LD with both causative mutations imparted a portion of the epistatic effect on these SNP, resulting in significant associations with the simulated traits.

In the case of genotyping, the ACA should select a subset of animals that, when genotyped, should give an optimal performance in terms of extrapolating the alleles of non-genotyped animals. Therefore, the objectives were to investigate the usefulness of a search algorithm as implemented by Ressom *et al.* (2006) to optimize the amount of information that can be ex‐ tracted from a pedigree while only genotyping a small portion. The results of the proposed method are compared to other viable methods to ascertain any potential gain. The proce‐ dures were tested using simulated pedigrees and actual beef cattle pedigrees of varying

*C.1 Ant colony optimization*: The ACA is initialized with all features having an equal baseline level of pheromone which is used to compute *Pm*(0) for all features. Using the PDF as de‐ fined in equation (1), each of *j* artificial ants will select a subset *Sk* of *n* features from the

Following the update of pheromone levels according to equation (2), the PDF is updated ac‐ cording to equation (1) and the process is repeated until some convergence criteria are met. Upon convergence the optimal subset of features is select based in the level of pheromone

In the specific case of selecting individuals for genotyping, the features are candidate ani‐ mals for genotyping from a full or partial pedigree. The pheromone of some feature, *m*, in the current study was proportional to the sum of an animal's number of mates and number

( ) *m m <sup>m</sup>*

*m m <sup>m</sup>*

Outside of actual ant colonies, and with regard in particular to the current study, it is diffi‐ cult to assign a biological explanation to the evaporation rate or*ρ*. Consequently, a relatively small value of 0.01 was chosen in an attempt to reach convergence faster. For each of *j* artifi‐ cial ants, a subset of animals was chosen equal to approximately 5% of the pedigree size.

For the five replicates of simulated pedigrees, 100 ants were used for each of 30,000 iter‐ ations. The evaporation rate was set equal to 0.01. The criterion used for evaluating can‐ didates was a function of their number of mates and number of offspring. Each animal in the pedigree was randomly assigned to be either homozygous or heterozygous. The probability of an animal being assigned to one of these two groups was dependent on

*t numoff nummate*

1

=

*n*

*m*

where *numoffm* and *nummatem* were the number of offspring and number of mates for animal *m* at time *t,* respectively. Consequently, the performance of a particular subset, *Sk*, is deter‐ mined the by the cumulative sum as described above for each of *n* animals in the subset.

*t numoff nummate* = + (9)

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

193

= + å (10)

t

( )

t

sizes and structures.

sample space *S* containing all features.

trail deposited on each feature.

of offspring

**Figure 4.** Plot of the cumulative distribution of SNP, identified as have significant associations when using a) ACA/G/D using 2 loci model (5.19%) b) RG/D, based on linkage disequilibrium with the causative mutations

### **2.3. Ant colony optimization as a method for strategic genotype sampling**

Interest in identifying QTL of economic importance for marker-assisted selection (MAS) in livestock populations has increased greatly in the past decade. Yet, it may not be via‐ ble to genotype each animal due to cost, time or lack of availability of DNA. A method that would allow for a selected sample (e.g. 5%) of the population to be genotyped and at the same time inferring with high probability genotypes for the remaining animals in the population could be beneficial. By using such a method, fewer animals in a popula‐ tion would be needed for genotyping which would decrease the time and cost of geno‐ typing. Theoretically the problem at hand is simple to solve. If it were possible to evaluate every possible subset of animals equal to the desired size (e.g. 5%) then the op‐ timal solution could be found. However, this is computationally impossible at the cur‐ rent time. Consequently a more feasible solution is needed. An intuitive solution would be one that selects animals based on their relationship with other animals in the pedi‐ gree. However, the heterozygosity and the structure of the pedigree play important roles as well. Consequently, the problem is one of optimization.

In the case of genotyping, the ACA should select a subset of animals that, when genotyped, should give an optimal performance in terms of extrapolating the alleles of non-genotyped animals. Therefore, the objectives were to investigate the usefulness of a search algorithm as implemented by Ressom *et al.* (2006) to optimize the amount of information that can be ex‐ tracted from a pedigree while only genotyping a small portion. The results of the proposed method are compared to other viable methods to ascertain any potential gain. The proce‐ dures were tested using simulated pedigrees and actual beef cattle pedigrees of varying sizes and structures.

To determine the effectiveness of the permutation on pheromone levels, the cumulative dis‐ tribution, based on LD with causative mutations, of SNP identified as being significantly as‐ sociated with simulated trait by ACA/G/D and RG/D were plotted and can be found in Fig. (4). Despite similarities in the average number of SNP identified by ACA/G/D (15.4) and RG/D (22), the distributions of these SNP, differed substantially. In contrast to RG/D, the ACA/G/D identified a large number of SNP having LD between.35-.45. These SNP corre‐ sponded to the broad plateau of SNP in LD with SNP 2041. Unlike RG/D, the ACA/G/D also identified several SNP (5.19%) having less than.10 LD with either of the causative mutations, an unexpected result given the strict family-wise significance thresholds (α=0.05) imposed on all models. Surprisingly, both methodologies identified a large number of SNP having LD of approximately ~.2. Upon closer examination it was found that these SNP had LD of ~. 2 with both causative mutations, likely artifacts of the data resulting from the relatively small sample size. The LD with both causative mutations imparted a portion of the epistatic

192 Ant Colony Optimization - Techniques and Applications

effect on these SNP, resulting in significant associations with the simulated traits.

(a) (b)

**Figure 4.** Plot of the cumulative distribution of SNP, identified as have significant associations when using a) ACA/G/D

Interest in identifying QTL of economic importance for marker-assisted selection (MAS) in livestock populations has increased greatly in the past decade. Yet, it may not be via‐ ble to genotype each animal due to cost, time or lack of availability of DNA. A method that would allow for a selected sample (e.g. 5%) of the population to be genotyped and at the same time inferring with high probability genotypes for the remaining animals in the population could be beneficial. By using such a method, fewer animals in a popula‐ tion would be needed for genotyping which would decrease the time and cost of geno‐ typing. Theoretically the problem at hand is simple to solve. If it were possible to evaluate every possible subset of animals equal to the desired size (e.g. 5%) then the op‐ timal solution could be found. However, this is computationally impossible at the cur‐ rent time. Consequently a more feasible solution is needed. An intuitive solution would be one that selects animals based on their relationship with other animals in the pedi‐ gree. However, the heterozygosity and the structure of the pedigree play important roles

using 2 loci model (5.19%) b) RG/D, based on linkage disequilibrium with the causative mutations

**2.3. Ant colony optimization as a method for strategic genotype sampling**

as well. Consequently, the problem is one of optimization.

*C.1 Ant colony optimization*: The ACA is initialized with all features having an equal baseline level of pheromone which is used to compute *Pm*(0) for all features. Using the PDF as de‐ fined in equation (1), each of *j* artificial ants will select a subset *Sk* of *n* features from the sample space *S* containing all features.

Following the update of pheromone levels according to equation (2), the PDF is updated ac‐ cording to equation (1) and the process is repeated until some convergence criteria are met. Upon convergence the optimal subset of features is select based in the level of pheromone trail deposited on each feature.

In the specific case of selecting individuals for genotyping, the features are candidate ani‐ mals for genotyping from a full or partial pedigree. The pheromone of some feature, *m*, in the current study was proportional to the sum of an animal's number of mates and number of offspring

$$
\pi\_m(t) = \text{numoff}\_m + \text{numumate}\_m \tag{9}
$$

where *numoffm* and *nummatem* were the number of offspring and number of mates for animal *m* at time *t,* respectively. Consequently, the performance of a particular subset, *Sk*, is deter‐ mined the by the cumulative sum as described above for each of *n* animals in the subset.

$$\pi\_m(t) = \sum\_{m=1}^{n} numoff\_{\mathbb{M}} + numrate\_m \tag{10}$$

Outside of actual ant colonies, and with regard in particular to the current study, it is diffi‐ cult to assign a biological explanation to the evaporation rate or*ρ*. Consequently, a relatively small value of 0.01 was chosen in an attempt to reach convergence faster. For each of *j* artifi‐ cial ants, a subset of animals was chosen equal to approximately 5% of the pedigree size.

For the five replicates of simulated pedigrees, 100 ants were used for each of 30,000 iter‐ ations. The evaporation rate was set equal to 0.01. The criterion used for evaluating can‐ didates was a function of their number of mates and number of offspring. Each animal in the pedigree was randomly assigned to be either homozygous or heterozygous. The probability of an animal being assigned to one of these two groups was dependent on the allelic frequencies such that if the allele frequencies were assumed to be 0.7/0.3 then approximately 58% of the animals would be categorized as homozygous based off of Hardy-Weinberg Laws of equilibrium. The assignment of homozygous/heterozygous sta‐ tus was performed each iteration. If a selected animal 2was homozygous then his/her number of mates and number of offspring were corrected such that for every homozy‐ gous offspring he/she had the number of offspring was corrected accordingly so that the number of offspring only reflected the number of heterozygous offspring. The same cor‐ rection was done for the number of mates. Similarly, if a selected animal was heterozy‐ gous, the number of offspring and the number of mates reflected a count of only homozygous individuals. An animal's probability of being selected was based off of maximizing the corrected sum of the animal's number of offspring and number of mates. The accuracy for evaluating a selected group of animals was proportional to this correct‐ ed sum. The uncorrected or original sum of each animal was used as prior information. Selected animals were chosen based off of their cumulative probability were assumed to have known genotypes for the peeling procedure. Simulated allele frequencies of 0.7/0.3 and 0.5/0.5 were used to assign genotypes to the animals in the pedigree.

where *n*1 and *n*2 were the number of animals with 2 and 1 allele(s) known and *na* was the total number of animals in the population. Furthermore, *n*1and *na* were multiplied by two

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

195

At the end of the peeling process those animals that had either one or two alleles known were retained for further analysis to determine the remaining unknown alleles in the popu‐ lation. In other words, those animals having one or two known alleles were used as prior information in the Gibbs sampling procedure for determining the remaining unknown al‐

*C.3 Gibbs sampling*: After the known alleles were determined by the peeling process descri‐ bed above, these alleles were used as prior information in the Gibbs Sampler to assign geno‐ types to the remaining animals in the population. For the base population animals, the unknown allele(s) were randomly sampled given the frequency of alleles in the population and the assumption of Hardy-Weinberg equilibrium. Unknown alleles for non-base popula‐ tion animals were randomly sampled from the parent's genotypes according to Mendelian rules. An equal weight was assumed for inheriting either the first or second allele from a parent. For a non-base population animal that had only one unknown allele, the unknown allele was sampled approximately half of the time from the sire's genotype and the remain‐ ing time from the dam's genotype. This was to compensate for incorrect assignment of the

At the end of the sampling process, a benefit function that described the total number of al‐ leles known in the population was computed. This function was computed from a combina‐ tion of known alleles and the probability of unknown alleles assigned during the sampling process. In order to be included in the benefit function, an allele in a particular position had to be equal to the true allele of the same position (i.e., *Bb* and *bB* were not equal). The proba‐

number of times ,

 was assigned ( ) . number of iterations *i j*

> 2 3 1 , ,1 ,2 1 1

*n n*

*i i Benefit n pa pa pa* = =

*a*

) and the number of known alleles, the benefit function was then computed as

2 [1 ( )] [ ( ) ( )],

) as previously defined. The percentage of alleles known after the Gibbs sampling

where*n*1, *n*2, and *n*3 were the number of animals with 2, 1 or 0 alleles known, respectively,

*i j i i*

, (*j* = 1 or 2) being assigned as the true allele j for animal *i* was calculated as:

*p a* = (12)

= ´+ + + + å å (13)

since each animal has two alleles.

known allele as illustrated in the above example.

,

*i j*

leles in the population.

bility of allele*ai*, *<sup>j</sup>*

Using *p*(*ai*, *<sup>j</sup>*

and *p*(*ai*, *<sup>j</sup>*

process, AKG, was such that

In the case of the real pedigree the same parameters were used as in the simulated pedigrees with the following exceptions; 100 ants were used for each of 5,000 iterations. The top 1,455 animals out of 29,101 were selected (5% of the total pedigree) based off of their cumulative probability were assumed to have known genotypes for the peeling procedure. In the case of the research beef cattle pedigree, 100 ants were used for each of 20,000 iterations. The top 434 out of 8,688 animals were selected (5% of the total pedigree) based on the same criteria.

*C.2 Peeling*: Given that genotypes in this study were assigned at random in the population, it is possible to extract additional genotypic information from the pedigree. Animals with missing genotypic information can be assigned one or both alleles given parental, progeny, or mate information. Given this trio of information sources and following an algorithm simi‐ lar to Qian and Beckmann (2002) and Tapadar et al. (2000), imputation on missing geno‐ types were made and additional genotypic information was garnered. For the current study it was assumed that there were no errors in the recorded pedigree resulting in all animals having known paternity and maternity. Whenever possible, maternal and paternal alleles were identified based on the inheritance. For the purpose of this study, the first allele was inherited from the sire and the second allele was inherited from the dam. If the parental ori‐ gin of an allele was unclear, then allele was arbitrarily assigned as either the paternal or ma‐ ternal allele.

After the peeling process, the number of animals with one or two alleles known was com‐ puted. This was done by simply counting the number of animals that were assigned either one or two alleles based on the peeling procedure described above. The percentage of alleles known based on the peeling procedure (AKP) was then computed as follows:

$$\mathbf{AK}\_{\mathrm{P}} = \left(\frac{(n\_1 \times 2) + n\_2}{n\_a \times 2}\right) \times 100,\tag{11}$$

where *n*1 and *n*2 were the number of animals with 2 and 1 allele(s) known and *na* was the total number of animals in the population. Furthermore, *n*1and *na* were multiplied by two since each animal has two alleles.

the allelic frequencies such that if the allele frequencies were assumed to be 0.7/0.3 then approximately 58% of the animals would be categorized as homozygous based off of Hardy-Weinberg Laws of equilibrium. The assignment of homozygous/heterozygous sta‐ tus was performed each iteration. If a selected animal 2was homozygous then his/her number of mates and number of offspring were corrected such that for every homozy‐ gous offspring he/she had the number of offspring was corrected accordingly so that the number of offspring only reflected the number of heterozygous offspring. The same cor‐ rection was done for the number of mates. Similarly, if a selected animal was heterozy‐ gous, the number of offspring and the number of mates reflected a count of only homozygous individuals. An animal's probability of being selected was based off of maximizing the corrected sum of the animal's number of offspring and number of mates. The accuracy for evaluating a selected group of animals was proportional to this correct‐ ed sum. The uncorrected or original sum of each animal was used as prior information. Selected animals were chosen based off of their cumulative probability were assumed to have known genotypes for the peeling procedure. Simulated allele frequencies of 0.7/0.3

194 Ant Colony Optimization - Techniques and Applications

and 0.5/0.5 were used to assign genotypes to the animals in the pedigree.

ternal allele.

In the case of the real pedigree the same parameters were used as in the simulated pedigrees with the following exceptions; 100 ants were used for each of 5,000 iterations. The top 1,455 animals out of 29,101 were selected (5% of the total pedigree) based off of their cumulative probability were assumed to have known genotypes for the peeling procedure. In the case of the research beef cattle pedigree, 100 ants were used for each of 20,000 iterations. The top 434 out of 8,688 animals were selected (5% of the total pedigree) based on the same criteria.

*C.2 Peeling*: Given that genotypes in this study were assigned at random in the population, it is possible to extract additional genotypic information from the pedigree. Animals with missing genotypic information can be assigned one or both alleles given parental, progeny, or mate information. Given this trio of information sources and following an algorithm simi‐ lar to Qian and Beckmann (2002) and Tapadar et al. (2000), imputation on missing geno‐ types were made and additional genotypic information was garnered. For the current study it was assumed that there were no errors in the recorded pedigree resulting in all animals having known paternity and maternity. Whenever possible, maternal and paternal alleles were identified based on the inheritance. For the purpose of this study, the first allele was inherited from the sire and the second allele was inherited from the dam. If the parental ori‐ gin of an allele was unclear, then allele was arbitrarily assigned as either the paternal or ma‐

After the peeling process, the number of animals with one or two alleles known was com‐ puted. This was done by simply counting the number of animals that were assigned either one or two alleles based on the peeling procedure described above. The percentage of alleles

1 2

(11)

( 2) AK 100, <sup>2</sup> *<sup>a</sup> n n n* æ ö ´ + = ´ ç ÷ ´ è ø

known based on the peeling procedure (AKP) was then computed as follows:

P

At the end of the peeling process those animals that had either one or two alleles known were retained for further analysis to determine the remaining unknown alleles in the popu‐ lation. In other words, those animals having one or two known alleles were used as prior information in the Gibbs sampling procedure for determining the remaining unknown al‐ leles in the population.

*C.3 Gibbs sampling*: After the known alleles were determined by the peeling process descri‐ bed above, these alleles were used as prior information in the Gibbs Sampler to assign geno‐ types to the remaining animals in the population. For the base population animals, the unknown allele(s) were randomly sampled given the frequency of alleles in the population and the assumption of Hardy-Weinberg equilibrium. Unknown alleles for non-base popula‐ tion animals were randomly sampled from the parent's genotypes according to Mendelian rules. An equal weight was assumed for inheriting either the first or second allele from a parent. For a non-base population animal that had only one unknown allele, the unknown allele was sampled approximately half of the time from the sire's genotype and the remain‐ ing time from the dam's genotype. This was to compensate for incorrect assignment of the known allele as illustrated in the above example.

At the end of the sampling process, a benefit function that described the total number of al‐ leles known in the population was computed. This function was computed from a combina‐ tion of known alleles and the probability of unknown alleles assigned during the sampling process. In order to be included in the benefit function, an allele in a particular position had to be equal to the true allele of the same position (i.e., *Bb* and *bB* were not equal). The proba‐ bility of allele*ai*, *<sup>j</sup>* , (*j* = 1 or 2) being assigned as the true allele j for animal *i* was calculated as:

$$p(a\_{i,j}) = \frac{\text{number of times } a\_{i,j} \text{ was assigned}}{\text{number of iterations}}.\tag{12}$$

Using *p*(*ai*, *<sup>j</sup>* ) and the number of known alleles, the benefit function was then computed as

$$Benefit = n\_1 \times 2 + \sum\_{i=1}^{n\_2} [1 + p(a\_{i,j})] + \sum\_{i=1}^{n\_3} [p(a\_{i,1}) + p(a\_{i,2})].\tag{13}$$

where*n*1, *n*2, and *n*3 were the number of animals with 2, 1 or 0 alleles known, respectively, and *p*(*ai*, *<sup>j</sup>* ) as previously defined. The percentage of alleles known after the Gibbs sampling process, AKG, was such that

$$\mathbf{AK}\_{\mathbf{G}} = \left(\frac{benefit}{n\_a \times \mathbf{2}}\right) \times 100 \,\mathrm{\,} \tag{14}$$

*C.5 Results of simulated pedigrees*: Table 9 presents results of the ACO and alternative meth‐ ods for analysis of the simulated pedigrees (Spangler 2008). The ant colony optimization method (ACO) appeared to be the most desirable method of those discussed in the current study. Compared to selecting 5% of the animals at random, ACO showed gains in AKP, AKG, and APTG ranging from 261.09 to 262.93%, 19.97 to 26.04%, and 23.5 to 29.6%, respec‐ tively. As compared to the favorable method of the alternative approaches, selecting males and females based of off the diagonal element of the inverse of the relationship matrix, the increase in AKP ranged from 4.98 to 5.16%. This gain is due to the amount of animals with both alleles known after the peeling process which was between 20.74 and 21.07% larger in favor of ACO. Admittedly, the gains in AKG were slight as compared to selecting males and females based of off the diagonal element of A-1, yet ACO still performed better. The in‐ crease in APTG ranged from 1.6 to 1.8% in favor of ACO over selecting males and females

Parameterb (0.30) (0.50) (0.30) (0.50) (0.30) (0.50) (0.30) ( 0.50)

2 alleles known 811.20 787.20 258.20 259.60 250.00 250.60 670.00 652.00 1 allele known 2,166.80 2,063.00 527.80 485.60 2,939.80 2,793.00 2,262.60 2,152.80 Benefit function 8,055.01 7,550.36 6,713.56 6,007.02 7,943.67 7,401.57 8,019.88 7,497.70 AKP 37.89 36.29 10.44 10.05 34.40 32.94 36.03 34.57 AKG 80.55 75.71 67.14 60.07 79.44 74.02 80.20 74.98 APTG 0.63 0.57 0.51 0.44 0.59 0.52 0.62 0.56

 Random= 5% selected at random, Males= 5% of males selected from their diagonal element of A-1, Males and fe‐ males= 2.5% males and 2.5% females selected from their diagonal element of A-1. Numbers in parenthesis are the true

**Table 9.** Number of animals with one or two alleles known, percentage of alleles known, and probability of assigning

*C.6 Real beef cattle pedigree*: Results from the ACO analysis can be found in Table 10 along with results from alternative approaches. The largest gains were seen in AKP which ranged from 150.00 to 171.62%, 2.95 to 3.04%, and from 1.80 to 1.94% as com‐ pared to random selection, selection of males and females from A-1, and selection of males from A-1, respectively. ACO also showed gains in AKG and APTG over random selection between 70.06 and 74.91% and between 14.3 and 15.4%, respectively. Table 3 shows advantages, although slight, of ACO over the methods using the diagonal element

*C.7 Research beef cattle pedigree*: Results from the ACO analysis and other approaches us‐ ing the same pedigree can be found in Table 11. As compared to randomly selecting 5%

allele frequencies used in the simulation. b Descriptions of the parameters can be found in equations 5-10

**ACO Random Males Males and females**

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

197

from their diagonal element.

the true genotype using other approachess

of A-1 for the parameters of AKG and APTG.

No. of animals with

a

where *benefit* was the benefit function computed above and *na* was the total number of ani‐ mals in the population.

During each round of the sampling process only one genotype of a given animal was as‐ signed as the true genotype. Thus, at the end of the sampling process every animal had a probability of having the true genotype,PTG*ig* , assigned as

$$\text{PTG}\_{\text{ig}} = \frac{\text{number of times genotype g was assigned}}{\text{total number of samples}},\tag{15}$$

where genotype *g* was the true genotype for animal*i*. The average probability of the true genotype being identified for every animal in the population (APTG) was computed using the following:

$$\text{APTG} = \frac{\sum\_{i=1}^{n\_a} \text{PTG}\_{ig}}{n\_a} \text{'} \tag{16}$$

where PTG*ig* was defined as above and *na* was the total number of animals in the popula‐ tion. In contrast to the benefit function, APTG only required that the animal have the correct genotype—*Bb* was considered the same genotype as *bB*—and therefore was able to compen‐ sate for the incorrect allele position and sampling the correct unknown allele.

*C.4 Simulation*: A simulation using an animal model was carried out to investigate two meth‐ ods of selecting animals for genotyping and two methods of maximizing the genetic infor‐ mation of the population. A pedigree with four over-lapping generations was simulated. The base population included 500 unrelated animals and subsequent generations consisted of 1,500 animals with a total of 5,000 animals generated. For the simulated pedigrees as well as the real pedigrees, one gene with two alleles was simulated for every animal in the pedi‐ gree file. Genotypes of the base population animals were assigned based on allele frequen‐ cies. For the subsequent generations, genotypes were randomly assigned using the parent's genotype, where an equal chance of passing either the first or second allele was assumed. Five replicates of the simulated data were generated.

Two different frequencies for the favorable allele were used in the simulation and analyses. The frequencies were 0.30, and 0.50. For the analyses using Gibbs sampling, a total chain length of 25,000 iterations of the Gibbs sampler was run, where the first 5,000 iterations were discarded as burn-in.

*C.5 Results of simulated pedigrees*: Table 9 presents results of the ACO and alternative meth‐ ods for analysis of the simulated pedigrees (Spangler 2008). The ant colony optimization method (ACO) appeared to be the most desirable method of those discussed in the current study. Compared to selecting 5% of the animals at random, ACO showed gains in AKP, AKG, and APTG ranging from 261.09 to 262.93%, 19.97 to 26.04%, and 23.5 to 29.6%, respec‐ tively. As compared to the favorable method of the alternative approaches, selecting males and females based of off the diagonal element of the inverse of the relationship matrix, the increase in AKP ranged from 4.98 to 5.16%. This gain is due to the amount of animals with both alleles known after the peeling process which was between 20.74 and 21.07% larger in favor of ACO. Admittedly, the gains in AKG were slight as compared to selecting males and females based of off the diagonal element of A-1, yet ACO still performed better. The in‐ crease in APTG ranged from 1.6 to 1.8% in favor of ACO over selecting males and females from their diagonal element.

AKG 100, <sup>2</sup> *<sup>a</sup> benefit n* æ ö = ´ ç ÷ ´ è ø

number of times genotype was assigned PTG , total number of samples *ig*

probability of having the true genotype,PTG*ig* , assigned as

mals in the population.

196 Ant Colony Optimization - Techniques and Applications

the following:

where *benefit* was the benefit function computed above and *na* was the total number of ani‐

During each round of the sampling process only one genotype of a given animal was as‐ signed as the true genotype. Thus, at the end of the sampling process every animal had a

where genotype *g* was the true genotype for animal*i*. The average probability of the true genotype being identified for every animal in the population (APTG) was computed using

> 1 PTG

> > *a n*

where PTG*ig* was defined as above and *na* was the total number of animals in the popula‐ tion. In contrast to the benefit function, APTG only required that the animal have the correct genotype—*Bb* was considered the same genotype as *bB*—and therefore was able to compen‐

*C.4 Simulation*: A simulation using an animal model was carried out to investigate two meth‐ ods of selecting animals for genotyping and two methods of maximizing the genetic infor‐ mation of the population. A pedigree with four over-lapping generations was simulated. The base population included 500 unrelated animals and subsequent generations consisted of 1,500 animals with a total of 5,000 animals generated. For the simulated pedigrees as well as the real pedigrees, one gene with two alleles was simulated for every animal in the pedi‐ gree file. Genotypes of the base population animals were assigned based on allele frequen‐ cies. For the subsequent generations, genotypes were randomly assigned using the parent's genotype, where an equal chance of passing either the first or second allele was assumed.

Two different frequencies for the favorable allele were used in the simulation and analyses. The frequencies were 0.30, and 0.50. For the analyses using Gibbs sampling, a total chain length of 25,000 iterations of the Gibbs sampler was run, where the first 5,000 iterations were

*ig*

APTG ,

*i*

<sup>=</sup> =

sate for the incorrect allele position and sampling the correct unknown allele.

Five replicates of the simulated data were generated.

discarded as burn-in.

*a n*

*<sup>g</sup>* <sup>=</sup> (15)

å (16)

(14)


a Random= 5% selected at random, Males= 5% of males selected from their diagonal element of A-1, Males and fe‐ males= 2.5% males and 2.5% females selected from their diagonal element of A-1. Numbers in parenthesis are the true allele frequencies used in the simulation. b Descriptions of the parameters can be found in equations 5-10

**Table 9.** Number of animals with one or two alleles known, percentage of alleles known, and probability of assigning the true genotype using other approachess

*C.6 Real beef cattle pedigree*: Results from the ACO analysis can be found in Table 10 along with results from alternative approaches. The largest gains were seen in AKP which ranged from 150.00 to 171.62%, 2.95 to 3.04%, and from 1.80 to 1.94% as com‐ pared to random selection, selection of males and females from A-1, and selection of males from A-1, respectively. ACO also showed gains in AKG and APTG over random selection between 70.06 and 74.91% and between 14.3 and 15.4%, respectively. Table 3 shows advantages, although slight, of ACO over the methods using the diagonal element of A-1 for the parameters of AKG and APTG.

*C.7 Research beef cattle pedigree*: Results from the ACO analysis and other approaches us‐ ing the same pedigree can be found in Table 11. As compared to randomly selecting 5% of the animals, ACO showed increases in AKP, AKG, and APTG ranging from 241.24 to 302.58%, 42.93 to 43.17%, and 20.9 to 38.0%, respectively. Realized gains in AKP of ACO over selecting males from A-1 or males and females from A-1 ranged from 8.78 to 10.15%, and 2.04 to 3.40%, respectfully.

**ACO Random Males Males and females**

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

199

Parameterb (0.30) (0.50) (0.30) (0.50) (0.30) (0.50) (0.30) ( 0.50)

2 alleles known 975.00 720.00 452.00 458.00 438.00 439.00 1,082.00 751.00 1 allele known 5,101.00 4,009.00 847.00 682.00 5,525.00 4,132.00 4,747.00 3,768.00 Benefit function 13,916.18 11,990.71 9,719.53 8,284.42 14,113.18 12,017.80 13,743.44 11,848.01 AKP 40.58 31.36 10.08 9.19 36.84 28.83 39.77 30.33 AKG 80.09 68.15 55.94 47.68 81.22 69.16 79.09 68.19 APTG 0.69 0.52 0.50 0.43 0.69 0.51 0.68 0.52

 Random= 5% selected at random, Males= 5% of males selected from their diagonal element of A-1, Males and fe‐ males= 2.5% males and 2.5% females selected from their diagonal element of A-1. Numbers in parenthesis are the true

Furthermore, pedigrees from field data or from research projects will also have innate struc‐ tural differences. Research projects may be limited by the size of the population and thus only use a small number of sires. In this scenario it would also be possible for higher rates of inbreeding and larger numbers of loops in a pedigree due to a large number of full sibs.

In the current study, the simulated pedigrees are composed of approximately 10% sires, while the large beef cattle pedigree and the small research beef cattle pedigree contain ap‐ proximately 16 and 7% sires, respectively. Intuitively, as the proportion of sires goes up, the number of offspring per sire goes down. This explains the similarity of the results between the simulated pedigrees and the small research pedigree. Thus, it is expected that the ACO algorithm will be far superior to other alternatives when very small (few hundred animals) pedigrees are considered or in situations where more than 5% of animals are genotyped due

Ant colony optimization offers a new and unique solution to the optimization problem of selecting individuals for genotyping. The heuristics used in the current study such as the number of ants, number of iterations, and the evaporation rate are unique only to the pedi‐ grees used in the current study. Each pedigree will offer a different structure and thus re‐

When applied to the high-dimensional data sets, the ant colony algorithm achieved higher prediction accuracies than all other feature selection methods examined. In contrast to previ‐ ous applications of optimization algorithms, the ant colony algorithm yielded high accura‐

allele frequencies used in the simulation. b Descriptions of the parameters can be found in equations 5-10.

**Table 11.** Number of animals with one or two alleles known, percentage of alleles known, and probability of

assigning the true genotype using other approaches from a real beef cattle research pedigreea

to reduction in animal with large diagonal elements in A-1.

quire a different set of parameters.

**3. Conclusions**

No. of animals with

a

The results suggest that ACO is the most desirable method of selecting candidates for genotyping, particularly after peeling (AKP). From these results it appears that the num‐ ber of offspring and the number of mates along with the homozygosity of the genotyped animals is critical in the selection process. Consequently, in application it will be critical to have good estimates of allele frequencies prior to implementing the genotype sam‐ pling strategy proposed in the current study. Differences in performance of ACO do ex‐ ist between the pedigrees explored in the current study. This is due to the proportion of sires and dams that have large numbers of offspring and/or mates. In the dairy industry, for example, there may be only a small number of sires in a pedigree but they may all be used heavily as in the case of the simulated pedigrees in the current study. In con‐ trast, a pedigree from the beef industry may have a larger proportion of sires but a large number of them may be used less frequently.


a Random= 5% selected at random, Males= 5% of males selected from their diagonal element of A-1, Males and fe‐ males= 2.5% males and 2.5% females selected from their diagonal element of A-1. Numbers in parenthesis are the true allele frequencies used in the simulation. b Descriptions of the parameters can be found in equations 5-10.

**Table 10.** Number of animals with one or two alleles known, percentage of alleles known, and probability of assigning the true genotype using other approaches from a real beef cattle pedigree a


a Random= 5% selected at random, Males= 5% of males selected from their diagonal element of A-1, Males and fe‐ males= 2.5% males and 2.5% females selected from their diagonal element of A-1. Numbers in parenthesis are the true allele frequencies used in the simulation. b Descriptions of the parameters can be found in equations 5-10.

**Table 11.** Number of animals with one or two alleles known, percentage of alleles known, and probability of assigning the true genotype using other approaches from a real beef cattle research pedigreea

Furthermore, pedigrees from field data or from research projects will also have innate struc‐ tural differences. Research projects may be limited by the size of the population and thus only use a small number of sires. In this scenario it would also be possible for higher rates of inbreeding and larger numbers of loops in a pedigree due to a large number of full sibs.

In the current study, the simulated pedigrees are composed of approximately 10% sires, while the large beef cattle pedigree and the small research beef cattle pedigree contain ap‐ proximately 16 and 7% sires, respectively. Intuitively, as the proportion of sires goes up, the number of offspring per sire goes down. This explains the similarity of the results between the simulated pedigrees and the small research pedigree. Thus, it is expected that the ACO algorithm will be far superior to other alternatives when very small (few hundred animals) pedigrees are considered or in situations where more than 5% of animals are genotyped due to reduction in animal with large diagonal elements in A-1.

Ant colony optimization offers a new and unique solution to the optimization problem of selecting individuals for genotyping. The heuristics used in the current study such as the number of ants, number of iterations, and the evaporation rate are unique only to the pedi‐ grees used in the current study. Each pedigree will offer a different structure and thus re‐ quire a different set of parameters.

## **3. Conclusions**

of the animals, ACO showed increases in AKP, AKG, and APTG ranging from 241.24 to 302.58%, 42.93 to 43.17%, and 20.9 to 38.0%, respectively. Realized gains in AKP of ACO over selecting males from A-1 or males and females from A-1 ranged from 8.78 to 10.15%,

The results suggest that ACO is the most desirable method of selecting candidates for genotyping, particularly after peeling (AKP). From these results it appears that the num‐ ber of offspring and the number of mates along with the homozygosity of the genotyped animals is critical in the selection process. Consequently, in application it will be critical to have good estimates of allele frequencies prior to implementing the genotype sam‐ pling strategy proposed in the current study. Differences in performance of ACO do ex‐ ist between the pedigrees explored in the current study. This is due to the proportion of sires and dams that have large numbers of offspring and/or mates. In the dairy industry, for example, there may be only a small number of sires in a pedigree but they may all be used heavily as in the case of the simulated pedigrees in the current study. In con‐ trast, a pedigree from the beef industry may have a larger proportion of sires but a large

Parameterb (0.30) (0.50) (0.30) (0.50) (0.30) (0.50) (0.30) ( 0.50)

2 alleles known 1,767.00 1,706.00 1,505.00 1,501.00 1,473.00 1,470.00 2,086.00 1,999.00

1 allele known 11,451.00 10,382.00 2,508.00 2,144.00 11,756.00 10,607.00 10,376.00 9,398.00

Benefit function 34,977.61 32,547.06 20,569.53 18,609.00 34,876.62 32,282.40 34,005.21 31,456.36

AKP 25.75 23.70 9.48 8.84 25.26 23.28 24.99 23.02

AKG 60.10 55.92 35.34 31.97 59.92 55.47 58.43 54.05

APTG 0.45 0.40 0.39 0.35 0.44 0.39 0.44 0.40

allele frequencies used in the simulation. b Descriptions of the parameters can be found in equations 5-10.

**Table 10.** Number of animals with one or two alleles known, percentage of alleles known, and probability of

assigning the true genotype using other approaches from a real beef cattle pedigree a

 Random= 5% selected at random, Males= 5% of males selected from their diagonal element of A-1, Males and fe‐ males= 2.5% males and 2.5% females selected from their diagonal element of A-1. Numbers in parenthesis are the true

**ACO Random Males Males and females**

and 2.04 to 3.40%, respectfully.

198 Ant Colony Optimization - Techniques and Applications

number of them may be used less frequently.

No. of animals with

a

When applied to the high-dimensional data sets, the ant colony algorithm achieved higher prediction accuracies than all other feature selection methods examined. In contrast to previ‐ ous applications of optimization algorithms, the ant colony algorithm yielded high accura‐ cies without the need to pre-select a small percentage of genes. Furthermore, the ant colony algorithm was able to identify small subsets of features with high predictive abilities and bi‐ ological relevance. In the presence of simulated epistasis, the proposed optimization meth‐ odology obtained substantial increases in power, demonstrating the effectiveness of machine learning approaches for the analysis of marker association studies in which gene interactions may be present. Although the ACA methods identified more SNP markers that could be construed as false positives, the use of a more stringent threshold eliminated the problem without greatly reducing the advantage of the ACA, in terms of power, when com‐ pared to other methods. The results of this study provide compelling evidence that the ACA is capable of efficiently modeling complex biological problems, such as the model proposed in this study.

[5] Barendse, W., Harrison, B. E., Hawken, R. J., Ferguson, D. M., Thompson, J. M., Tho‐ mas, M. B., and R. J. Bunch. 2007. Epistasis between Calpain 1 and its inhibitor Cal‐

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

201

[6] Coutinho, A. M., Sousa, I., Martins, M. et al. 2007. Evidence for epistasis between SLC6A4 and ITGB3 in autism etiology and in the determination of platelet serotonin

[7] Ding, Y. P., Wu, Q. S., and Q. D. Su. 2005. Multivariate Calibration Analysis for metal porphyrin mixtures by an ant colony algorithm. Analytical Sciences. 21:327-330.

[8] Dorigo M., Di Caro G. & Gambardella L.M. (1999) Ant algorithms for discrete opti‐

[9] Dorigo, M. and L. M. Gambardella. 1997. Ant colonies for the travailing salesman

[10] Golub,T.R., Slonim,D.K., Tomayo,P., Huard,C., Gaasenbeek,M., Mesirov,J.P., Col‐ ler,H., Loh,M.L., Downing,J.R., Caligiuri,M.A., Bloomfield,C.D. and E. S. Lander (1999) Molecular classification of cancer: class discovery and class prediction by gene

[11] Gonzalez, J. R., Armengol, L., Sole, X., Guino, E., Mercader, J. M., Estivill, X., and V. Moreno. 2007. SNPassoc: an R package to perform whole genome association studies.

[12] Hong,J. and S. Cho (2006) Efficient huge-scale feature with speciated genetic algo‐

[13] Hugot, J. P., Chamaillard, M., Zouali, H. et al. 2001. Association of NOD2 leucinerich repeat variants with susceptibility to Crohn's disease. Nature. 411:599-603.

[14] Jefferey,I.B., Higgins,D.G. and A. Culhane (2006) Comparison and evaluation of methods for generating differentially expressed gene lists from microarray data,

[15] Lin,T., Liu,R., Chen,C., Choa,Y. and S. Chen (2006) Pattern classificationin DNA mi‐

[16] Liu,J.J., Cutler,G., Li,W., Pan,Z., Peng,S., Hoey,T., Chen,L. and X. B. Ling (2005) Mul‐ ticlass cancer classification and biomarker discover using GA-based algorithms *Bioin‐*

[17] Marchini, J., Donnelly, P., and L. R. Cardon. 2005. Genome-wide stregies for detect‐ ing multiple loci that influence complex diseases. Nat. Genetics. 37:413-417.

[18] Nagata,K., Kawajiri,A., Matsui,S., Takagishi,M., Shiromizu,T., Saitoh,N., Izawa,I. Kiyono,T., Itoh,T.J., Hotani,H. and M. Inagaki (2003) Filament formation of MSF-A, a mammalan Septin, in human mammary epithelial cells depends on interactions with

croarray data of multiple tumor types *Pattern Recognition,* 39, 2426-2438.

pastatin within breeds of cattle. Genetics 176:2601-2610.

levels. Hum. Genet. 121:243-256.

mization. Artificial Life 5, 137–72.

problem. BioSystems. 43:73-81.

Bioinformatics. 23(5):644-645

*BMC Bioinformatics*, 7.

*formatics,* 21, 2691-2697.

expression monitoring *Science,* 286, 531-537.

rithm *Pattern Recognition Lett*., 27, 143-150.

microtubules *J. of Biol. Chem.,* 278, 18538-18543

## **Author details**


## **References**


[5] Barendse, W., Harrison, B. E., Hawken, R. J., Ferguson, D. M., Thompson, J. M., Tho‐ mas, M. B., and R. J. Bunch. 2007. Epistasis between Calpain 1 and its inhibitor Cal‐ pastatin within breeds of cattle. Genetics 176:2601-2610.

cies without the need to pre-select a small percentage of genes. Furthermore, the ant colony algorithm was able to identify small subsets of features with high predictive abilities and bi‐ ological relevance. In the presence of simulated epistasis, the proposed optimization meth‐ odology obtained substantial increases in power, demonstrating the effectiveness of machine learning approaches for the analysis of marker association studies in which gene interactions may be present. Although the ACA methods identified more SNP markers that could be construed as false positives, the use of a more stringent threshold eliminated the problem without greatly reducing the advantage of the ACA, in terms of power, when com‐ pared to other methods. The results of this study provide compelling evidence that the ACA is capable of efficiently modeling complex biological problems, such as the model proposed

, S. Smith1

1 Department of Animal and Dairy Science, The University of Georgia, Athens, Greece

[1] Akiyama,T. and Y Kawasaki (2006) Wnt signaling and the actin cytoskeleton *Onco‐*

[2] Albrecht, A., Vinterbo,S.A. and L. O. Machado 2003, 'An epicurean learning ap‐ proach to gene-expression data classification', *Artif. Intell in Medicine*, 28, 75-87.

[3] Antonov,A.V., Tetko,I.V., Mader,M.T., Budczies,J. and H. W. Mewes (2004) Optimi‐ zation models for cancer classification: extracting gene interaction information from

[4] Bagirov,A.M., Ferguson,B., Ivkovic,S., Saunders,G. and J. Yearwood (2003) New al‐ gorithms for multi-class cancer diagnosis using tumor gene expression signatures *Bi‐*

, E. H. Hay1

and K. Bertrand1

in this study.

**Author details**

**References**

R. Rekaya1,2,3\*, K. Robbins4

200 Ant Colony Optimization - Techniques and Applications

, M. Spangler5

2 Department of Statistics, The University of Georgia, Athens, Greece

3 Institute of Bioinformatics, The University of Georgia, Athens, Greece

5 Animal Science Department, University of Nebraska, Lincoln, NE, USA

microarray expression data *Bioinformatics,* 20, 644-652.

\*Address all correspondence to: rrekaya@uga.edu

4 Dow AgroSciences, Indianapolis, IN, USA

*gene,* 25, 7538-7544.

*oinformatics,* 19, 1800-1807.


[19] Ooi,C.H. and P. Tan (2003) Genetic algorithms applied to multi-class prediction for the analysis of gene expression data *Bioinformatics,* 19, 37-44.

[32] Spangler, M. L., K. R. Robbins, J. K. Bertrand, M. MacNeil, and R. Rekaya. 2008. Ant colony optimization as a method for strategic genotype sampling. Animal Genetics

Ant Colony Algorithm with Applications in the Field of Genomics

http://dx.doi.org/10.5772/52051

203

[33] Subramani,P., Sahu,R. and S. Verma, Feature selection using Haar wavelet power

[34] Tapadar P., Ghosh S. & Majumder P.P. (2000) Haplotyping in pedigrees via a genetic

[35] West M. (2003) Bayesian factor regression models in the "Large p, Small n" paradigm,

[36] Woon , P. Y., Kaisaki, P. J., Braganca, J., Bihoreau, M. T., Levy, J. C., Farrall, M., and D. Gauguir. 2007. Aryl hydrocarbon receptor nuclear translocator-like (BMAL1) is as‐ sociated with susceptibility to hypertension and type 2 diabetes. Proc. Natl. Acad.

40: 308 – 314.

spectrum *BMC Bioinformatics,* 7:432.

algorithm. Human Heredity 50, 43–56.

*Bayesian Statistics*, 7, 723-732.

Sci. 104(36):14412-14417.


[32] Spangler, M. L., K. R. Robbins, J. K. Bertrand, M. MacNeil, and R. Rekaya. 2008. Ant colony optimization as a method for strategic genotype sampling. Animal Genetics 40: 308 – 314.

[19] Ooi,C.H. and P. Tan (2003) Genetic algorithms applied to multi-class prediction for

[20] Peng,S., Xu,Q., Ling,X.B., Peng,X., Du,W. and L. Chen Molecular classification of can‐ cer types from microarray data using the combination of genetic algorithms and sup‐

[21] Qian D. & Beckmann L. (2002) Minimum-recombinant haplotyping in pedigrees.

[22] Ramaswamy,S., Tamayo,P., Rifkin,R., Mukherjee,S., Yeang,C., Angelo,M., Ladd,C., Reich,M., Latulippe,E., Mesirov,J.P., Poggio,T., Gerald,W., Loda,M., Lander,E.S. and T. R. Golub (2001) Multiclass cancer diagnosis using tumor gene expression signa‐

[23] Rekaya, R, K. Robbins. (2009). Ant colony algorithm for analysis of gene interaction in high-dimensional association data. Revista Brasileira de Zootecnia. doi: 10.1590/

[24] Ressom,H.W., Varghese,R.S., Orvisky,E., Drake,S.K., Hortin,G.L., Abdel-Hamid,M. Loffredo,C.A. and R. Goldman (2006) Ant colony optimization for biomarker identi‐ fication from MALDI-TOF mass spectra *Proc. of the 28th EMBS Annual Inter. Conf.,*

[25] Robbins, K. R., Zhang, W., R. Rekaya, and J. K. Bertrand. 2007. The use of the ant col‐ ony algorithm for analysis of high-dimension gene expression data sets. 58th Annual

[26] Robbins, K. R., Zhang, W., J. K. Bertrand, R. Rekaya. 2008. Ant colony optimization for feature selection in high dimensionality data sets. Math Med Biol. 24(4):413-426.

[27] Robbins K, K. Bertrand, and R. Rekaya. 2011. The use of the ant colony algorithm for the detection of marker associations in the presence of gene interactions. Internation‐

[28] Scott,M., McCluggage,W.G., Hillan,K.J., Hall,P.A. and S. E. H. Russell (2006) Altered patterns of transcription of th septin gene, SEPT9, in ovarian tumorgenesis *Int. J. Can‐*

[29] Shen,R., Ghosh,D., Chinnaiyan,A. and Z. Meng Eigengene-based linear discriminant model for tumor classification using gene expression microarray data *Bioinformatics,*

[30] Shymygelska, A. and H. H. Hoos. 2005. An ant colony optimization algorithm for the 2D and 3D hydrocarbon polar protein folding program. BMC Bioinformatics. 6:30.

[31] Sinnwell, J. P. and D. J. Schaid. 2005. haplo.stats: Statistical Analysis of Haplotypes with Traits and Covariates when Linkage Phase is Ambiguous. R package version

Meeting of the European Association for Animal Production (EAAP):167.

the analysis of gene expression data *Bioinformatics,* 19, 37-44.

port vector machines *FEBS Letters,* 555, 358-362.

American Journal of Human Genetics 70, 1434–45.

al Journal of Bioinformatics Research, 2:227-235.

tures *PNAS,* 98, 15149-15154.

202 Ant Colony Optimization - Techniques and Applications

S1516-35982009001300011.

4560-4563.

*cer,* 118, 1325-1329.

22, 2635-2642.

1.2.2.


## *Edited by Helio J.C. Barbosa*

Ant Colony Optimization (ACO) is the best example of how studies aimed at understanding and modeling the behavior of ants and other social insects can provide inspiration for the development of computational algorithms for the solution of difficult mathematical problems. Introduced by Marco Dorigo in his PhD thesis (1992) and initially applied to the travelling salesman problem, the ACO field has experienced a tremendous growth, standing today as an important nature-inspired stochastic metaheuristic for hard optimization problems.

This book presents state-of-the-art ACO methods and is divided into two parts: (I) Techniques, which includes parallel implementations, and (II) Applications, where recent contributions of ACO to diverse fields, such as traffic congestion and control, structural optimization, manufacturing, and genomics are presented.

Ant Colony Optimization - Techniques and Applications

Ant Colony Optimization

Techniques and Applications

*Edited by Helio J.C. Barbosa*

Photo by Mirexon / iStock