Abstract
This paper deals with usage of an alternative tool for symbolic regression—analytic programming which is able to solve various problems from the symbolic domain, as well as genetic programming and grammatical evolution. This paper describes a setting of an optimal trajectory for a robot (originally designed as an artificial ant on Santa Fe trail) solved by means of analytic programming. Firstly, main principles of analytic programming are described and explained. The second part shows how analytic programming was used for the application of finding a suitable trajectory step by step. Because analytic programming needs evolutionary algorithms for its run, three evolutionary algorithms were used—self-organizing migrating algorithm, differential evolution, and simulated annealing—to show that anyone can be used. The total number of simulations was 150 and results show that the first two used algorithms were more successful than not so robust simulated annealing.
1. Introduction
The term “symbolic regression” represents a process during which measured data is fitted and a suitable mathematical formula is obtained in an analytical way. This process is well known for mathematicians. It is used when a mathematical model of unknown data is needed. For long time, symbolic regression was a domain of humans but in the last few decades, computers have gone to foreground of interest in this field. Firstly, the idea of symbolic regression done by means of computer was proposed by Koza in genetic programming (GP) [1–3]. The other two approaches are grammatical evolution (GE) developed by Ryan et al. [4–6] and here described analytic programming (AP) designed in [7–9].
Genetic programming was the first tool for symbolic synthesis of the so-called programs done by means of computer instead of humans. The main idea comes from genetic algorithms (GAs) [10], which Koza uses in his GP. The ability to solve very difficult problems was proved many times, and hence, GP today can be applied, for example, to synthesize highly sophisticated electronic circuits, robot trajectory, biochemistry problems, and many others [2].
The other tool is GE which was developed in the last decade of 20th century by Conor Ryan. Gramatical evolution has one advantage compared to GP which is the ability to use arbitrary programming language not only LISP as in the case of the cannonical version of GP. In contrast to other evolutionary algorithms, GE was used with a few search strategies with a binary representation of the populations [5], as well as with other algorithms like those in [11, 12]. Other 2 interesting investigations using symbolic regression were carried out by Johnson [13] working on artificial immune systems and probabilistic incremental program evolution (PIPE), the work in [14] generates functional programs from an adaptive probability distribution over all possible programs.
This contribution demonstrates the use of a method which is independent of computer platform, programming language, and can use any evolutionary algorithm (as demonstrated in [7–9]) to find an optimal solution of the required task.
2. Analytic Programming
2.1. Description
Basic principles of the AP were developed in 2001. Until that time, mainly GP and GE existed. GP uses genetic algorithms while AP can be used with any evolutionary algorithm, independently of individual representation. To avoid any confusion based on the use of names according to the used algorithm, the name Analytic programming was chosen, because AP stands for synthesis of analytical solution by means of evolutionary algorithms [7–9].
AP was inspired, in general, by numerical methods in Hilbert spaces and by GP. Principles of AP [9] are somewhere between these two philosophies. From GP, an idea of evolutionary creation of symbolic solutions is taken into AP while from Hilbert spaces, an idea of synthesis of more complicated functions from elementary functions is adopted into AP. Analytic programming as well as GPis based on the set of functions, operators, and so-called terminals, which are usually constants or independent variables like
(i)functions: sin, tan, And, Or, and so forth,(ii)operators: +, −, *, /, dt, and so forth,(iii)terminals: 2.73, 3.14, , and so forth.
All these “mathematical” objects create a set which AP tries to synthesize the appropriate solution from. The set of mathematical objects are functions, operators, and so-called terminals (usually constants or independent variables). All these objects are mixed together as shown in Figure 1and consist of functions with different number of arguments. Because of the variability of the content of this set, it is called for article purposes general functional set (GFS). The structure of GFS is nested, that is, it is created by subsets of functions according to the number of their arguments. The content of GFS is dependent only on the user. Various functions and terminals can be mixed together. For example, is a set of all functions, operators, and terminals, is a subset containing functions with only three arguments, represents only terminals, and so forth.
This nested structure is necessary that the main principle of AP can work without any difficulties. The core of AP is based on discrete set handling, proposed in [15, 16] (see Figure 2). Discrete set handling (DSH) shows itself as a universal interface between EA and symbolically solved problem. That is why AP can be used almost by any evolutionary algorithm.
Briefly said, DSH works with integer indexes which represent numerical or nonnumerical expressions (operators, functions, etc.) in a discrete set. This index then serves like a pointer into a discrete set. Based on that, appropriate objects are chosen for cost function evaluation [16]. During an evolutionary process, only indexes are used for all evolutionary operations. Objects from the discrete set are used (by means of integer index) only in cost function, whereas according to integer index, a symbolic structure is synthesized and consequently evaluated.
2.2. Mapping Method in AP
The nested structure presence in GFS is vitally important for AP. It is used to avoid synthesis of pathological programs, that is, programs containing functions without arguments, and so forth. Performance of AP is, of course, improved if functions of GFS are expertly chosen based on experiencies with solved problem.
The important part of the AP is a sequence of mathematical operations which are used for program synthesis. These operations are used to transform an individual of a population into a suitable program. Mathematically said, it is mapping from an individual domain into a program domain. This mapping consists of two main parts. The first part is called discrete set handling (DSH) and the second one is security procedures which do not allow synthesizing pathological programs.
Discrete set handling proposed in [15, 16] is used to create an integer index, which is used in the evolutionary process like an alternative individual handled in EA by method of integer handling. The method of DSH, when used, allows handling arbitrary objects including nonnumeric objects like linguistic terms (hot, cold, dark, etc.), logic terms (true, false), or other user defined functions. In the AP, DSH is used to map an individual into GFS and together with security procedures (SP) creates the aforementioned mapping which transforms arbitrary individual into a program. Individuals in the population consist of integer parameters, that is, an individual is an integer index pointing into GFS.
Analytic programming is basically a series of function mapping. Figure 3 demonstrated an artificial example of how a final function is created from an integer individual. Number 1 in the position of the first parameter of integer index means that the operator “+” from is used. Because the operator “+” has to have at least two arguments, next two index pointers 6 (sin from ) and 7 (cos from ) are dedicated to this operator as its arguments. Both functions, sin and cos, are one-argument functions so the next unused pointers 8 (tan from ) and 9 ( from ) are dedicated to sin and cos function. Because as an argument of cos variable is used, this part of resulting function is closed ( is zero-argument) in its AP development. One-argument function tan remains, and because there is one unused pointer 9 tan is mapped on “” which is on the 9th position in GFS.
To avoid synthesis of pathological functions a few security “tricks” are used in AP. The first one is that GFS consists of subsets containing functions with the same number of arguments. Existence of this nested structure is used in the special security subroutine which is measuring how far the end of individual is, and according to it, objects from different subsets are selected to avoid pathological function synthesis. Precisely, if more arguments are desired than possible (the end of the individual is near), function will be replaced by other function with the same index pointer from subset with lower number of arguments. For example, it may happen that the last argument for one argument function will not be a terminal (zero-argument function). If pointer is bigger than length of subset, that is, the pointer is 5 and is used , then the element is selected according to element = pointer_value mod number_of_elements_in_. In this example, case-selected element would be variable t (see in Figure 1).
GFS needs to be constructed not only from clear mathematical functions as demonstrated but also from other user-defined functions, which can be used, for example, logical functions, functions which represent elements of electrical circuits, or robot movement commands.
2.3. Versions of AP
Today, AP exists in three versions: , , and . In all three versions, the same sets of functions, terminals, and so forth, as Koza use in GP [1–3] are necessary for the program synthesis. works as described earlier and the formulas do not contain any constants. The second version () is modified in the sense of constant estimation. For example, when Koza uses randomly generated constants in the so-called sextic problem [3], AP uses only one (), which is inserted into the formula at various places by evolutionary processing. The function can look as follows: When the program is synthesized, then all “” are indexed so that , are obtained from (2), and then all are estimated by second evolutionary algorithm, and the result is in(3): Because EA “works under” EA (i.e., program estimation of ), this version is called AP with metaevolution–-. As this version was quite time-consuming, another modification of was done extending the second version by estimation of . It is done by suitable methods of nonlinear fitting (). This method has shown the most promising performance when unknown constants are present.
2.4. Security Procedures
Security procedures (SPs) are in the AP as well as in GP, used to avoid various critical situations. In the case that AP security procedures were not developed for AP purposes after all, but they are mostly integrated parts of AP. However sometimes they have to be defined as a part of cost function, based on kind of situation (e.g., situation 2, 3, and 4, etc., see what follows). Critical situations are like
(1)pathological function (e.g., without arguments, self-looped),(2)functions with imaginary or real part (if not expected),(3)infinity in functions (e.g., dividing by 0),(4)“frozen” functions (e.g., extremely long time to get a cost value: hours).
Simply as an SP can be regarded here mapping from an integer individual to the program which is checked for how far the end of the individual is, and based on this information, a sequence of mapping is redirected into a subset with lower number of arguments. This satisfies that no pathological function will be generated. Another activities of SP are integrated part of cost function to satisfy items 2–4, and so forth.
2.5. Similarities and Differences
Because AP was partly inspired by GP, then between AP, GP, and GE are some differences as well as some logical similarities. A few of the most important ones are as follows.
I. Similarity
(i)Synthesized programs: AP as well as G0P and GE is able to do symbolic regression in general point of view. It means that output of AP is according to all important simulations [7–9] similar to programs from GP and GE (see http://www.fai.utb.cz/people/zelinka/ap).(ii)Functional set: operates in principle on the same set of terminals and functions as GP or GE.
II. Differences
(i) or use universal constant (difference) which is indexed after program synthesis.(ii)Individual coding: coding of an individual is different. Analytic programming uses an integer index instead of direct representation as in canonical GP. Grammatical evolution uses binary representation of an individual, which is consequently converted into integers for mapping into programs by means of BNF [4].(iii)Individual mapping: AP uses discrete set handling, [13] while GP in its fundamental form uses direct representation in Lisp [1] and GE uses grammar-Backus-Naur form (BNF) [4].(iv)Constant handling: GP uses a randomly generated subset of numbers, constants, GE utilises user-determined constants and AP uses only one constant for and , which is estimated by other EA or by nonlinear fitting.(v)Security procedures: to guarantee synthesis of nonpathological functions, procedures are used in AP which redirect the flow of mapping into subsets of a whole set of functions and terminals according to the distance to the end of the individual. If a pathological function is synthesized in GP, then synthesis is repeated. In the case of GE, when the end of an individual is reached, then mapping continues from the individual beginning, which is not the case of AP. It is designed so that a nonpathological program is synthesized before the end of the individual is reached (maximally when the end is reached).
2.6. Selected Solved Problems
During AP development and research simulations, a lot of various kinds of programs have been synthesized. In (2) a mathematical formula is shown to demonstrate complexity of synthesized formulas, which were randomly generated amongst 1000 formulas to check if the final structure is free of pathologies (i.e., if all functions have the right number of arguments, etc.). In this case, no attention was paid to mathematical reasonability of the following test programs based on clear mathematical functions. In what mentioned earlier, a different approach to the symbolic regression called analytic programming was described. Based on its results and structure, it can be stated that AP seems to be a universal candidate for symbolic regression by means of different search strategies. Problems on which AP was utilised were selected from test and theory problems domain as well as from real-life problems and are shown in following examples.
(i)Random synthesis of function from GFS, 1000 times repeated: the aim of this simulation was to check if pathological function can be generated by AP. In this simulation, randomly generated individuals were created and consequently transformed into programs and checked for their internal structure. No pathological program was identified [7].(ii)sin() approximation was repeated 100 times. Here AP was used to synthesize the program function sin() fitting [7].(iii) approximation was repeated 100 times, the same as in the previous example. Main aim was again fitting of dataset generated by a given formula [7].(iv)Solving of ordinary differential equations (ODE): ’’() = cos(), (0) = 1, , ’(0) = 0, ’() = 0, was repeated 100 times, in that case AP was looking for suitable function, which would solve this case of ODE [7].(v)Solving of ODE: (4 + )’’()’’ + 600() = 5000(), (0) = 0, (1) = 0, ’’(0) = 0, ’’(1) = 0, was repeated 5 times (due to longer time of simulation in the Mathematica environment). Again as in the previous case, AP was used to synthesize a suitable function-solution of this kind of ODE. This ODE was used from and represents a civil-engineering problem in reality [7].(vi)Boolean even and symmetry problems according to [1] for comparative reasons [9].(vii)Sextic and Quintic problems [8].(viii)Simple neural network synthesis by means of AP: a simple few layered NN synthesis was tested by AP [17].
Such elementary objects are usually simple mathematical operators (+, −, *, …), simple functions (sin, cos, And, Nor, etc.), user-defined functions, and so forth. Output of symbolic regression is a more complex “object” (formula, function, command, etc.), solving a given problem like data fitting of so-called Sextic and Quintic problem described by (4), [2, 8], randomly synthesized function (5) [8], as well as Boolean problems of parity and symmetry solution (basically logical circuits synthesis) (6) [2, 9]. However, (4)–(6) mentioned here are just only a few samples of numerous repeated experiments done by AP and are used to demonstrate how complex structures can be produced by symbolic regression in general sense for different problems: The rest of this article is an investigation on evolutionary synthesis of robot commands, which is well known in genetic programming as a Santa Fe trail for an artificial ant.
3. Problem Design
3.1. Santa Fe Description
The Santa Fe trail, demonstrated in Figure 4, was chosen from [18] to make a comparative study with the same problem which was solved by Koza in genetic programming [1].
The aim of the task is that an artificial ant should go through defined trail and eat all food which is there. From a simple point of view, it can be looked at it as on robot movements on some trail. Robot trajectory is, of course, very complex task but the more complex behaviour can be added later in further simulations.
The Santa Fe trail is defined as a field where food is set out. In Figure 4, a black field is food for the ant. The gray one is basically the same as a white field but, for clarity, was used the gray color. The gray fields represent obstacles (fields without food on the road) for the ant. If there would not be these holes, the ant could go directly through the way. It would be enough to go and see before ant if there is food. If yes, ant would go straight and eat the bait. If not, it would turn around and see where food is, and the cycle would repeat till the ant would eat the last bait.
In the real world, robots have obstacles in their moving. Therefore, also in this case, such approach was chosen. The first problem which ant has to overcome is the simple hole (position (8,27) in Figure 4). Second one is the two holes in the line (positions (13,16) and (13,17), or three holes (17,15), (17,16), (17,17). Next problem is the holes in the corners: one (position (13,8), two (1,8), (2,8), and three holes (17,15), (17,16), (17,17).
3.2. Set of Functions
The set of functions used for movements of the ant is as follows. As a set of variables , that is, in the case of this article there are functions which provide movements of an ant, without any argument which could be add during the process of evolution.
The set consist of
(i) = Left, Right, Move, where
: a set of variables and terminals, zero argument functions ,
Left: function for turning around in the anticlockwise direction,
Right: function for turning around in the clockwise direction,
Move: function for moving straight and if bait is in the field where the ant is moved, it is eaten.
This set of functions is not enough to make successfully a desired task. More functions are necessary, then a and a were set up:
(ii) = IfFoodAhead, Prog2,(iii) = Prog3, where the number in GFS means the arity of the functions inside, that is, the number of arguments which are needed to be evaluated correctly. Arguments are added to those functions during evolution process, as mentioned earlier in the description of AP.
IfFoodAhead is a decision function: the ant controls the field in front of it, and if there is food, the function in the field for truth argument is executed; otherwise, function in false position is performed.
Prog2 and Prog3 are the same function in the principle. They do 2 or 3 functions in the same time. These two functions were originally defined also in Koza’s approach but in AP, it is necessary because of the structure of generating the program.
3.3. Fitness Function
The aim of the ant is to eat all food on the way. There are 89 baits. This is so called raw fitness, and the value of cost function (7) is calculated as a difference between raw fitness and a number of baits eaten by an ant [1], which went through the grid according to just generated way: where Number_of_Food is number of eaten baits by an ant according to synthesized way.
The aim is to find such formula whose cost value is equal to zero. To obtain an appropriate solution, two constraints should be set up into a cost function. One is a limitation concerned to the number of steps. It is not desired to the ant to go field by field in the grid. A requirement to the fastest and the most effective way is desired. Then a limit of steps was equal to 600. According to the original assignment, 400 steps should be sufficient, but as the work in [19], Koza’s optimal solution was as in (8). However, as simple solution showed, 545 steps are necessary for an ant to eat all food in the Santa Fe trail. Functionality of (8) can be described in follows. If bait is in front of the ant, it moves on the field and eats the food. If there is nothing, it does the following 3 commands. If food is in front of the ant it moves and eats the food, if not it turns twice right. Next Prog2(Left, Right) is not necessary there, this is the reason why all program takes 545 steps instead of 404 in the case of no Prog2(Left, Right). Then next control of food in front of ant is again, if yes ant moves and eats the food. If not it turns left to the original direction as it was at the beginning of the program. If the cycle is somewhere interrupt (e.g., in the case of truth in the first function IfFoodAhead), the cycle is repeated still from the beginning until all food is not eaten or constrained steps are not reached.
The second constraint could be concerned to the length of the list of commands for an ant. The longer can cause the more steps to reach all food is eaten. In this preliminary study, this constraint was not set up, but in further studies, a penalization concerned this constraint will be surely used.
4. Used Evolutionary Algorithm
In this paper, self-organizing migrating algorithm (SOMA), differential evolution (DE), and simulated annealing (SA) were used as an evolutionary algorithm. For detailed information, see [15, 20, 21].
4.1. Differential Evolution (DE)
Differential evolution is a population-based optimization method that works on real-number-coded individuals [20]. For each individual in the current generation G, DE generates a new trial individual by adding the weighted difference between two randomly selected individuals and to a third randomly selected individual . The resulting individual is crossed-over with the original individual . The fitness of the resulting individual, referred to as a perturbated vector , is then compared to the fitness of . If the fitness of is greater than the fitness of , is replaced with , otherwise remains in the population as .
Differential evolution is robust, fast, and effective with global optimization ability. It does not require that the objective function is differentiable, and it works with noisy, epistatic, and time-dependent objective functions.
4.2. Self-Organizing Migrating Algorithm (SOMA)
SOMA is a stochastic optimization algorithm that is modeled on the social behaviour of cooperating individuals [15] It was chosen because it has been proven that the algorithm has the ability to converge towards the global optimum [15]. SOMA works on a population of candidate solutions in loops called migration loops. The population is initialized randomly distributed over the search space at the beginning of the search. In each loop, the population is evaluated and the solution with the highest fitness becomes the leader L. Apart from the leader, in one migration loop, all individuals will traverse the input space in the direction of the leader. Mutation, the random perturbation of individuals, is an important operation for evolutionary strategies (ESs). It ensures the diversity among the individuals, and it also provides the means to restore lost information in a population. Mutation is different in SOMA compared with other ES strategies. SOMA uses a parameter called PRT to achieve perturbation. This parameter has the same effect for SOMA as mutation has for GA.
The novelty of this approach is that the PRT Vector is created before an individual starts its journey over the search space. The PRT Vector defines the final movement of an active individual in search space.
The randomly generated binary perturbation vector controls the allowed dimensions for an individual. If an element of the perturbation vector is set to zero, then the individual is not allowed to change its position in the corresponding dimension.
An individual will travel a certain distance (called the PathLength) towards the leader in n steps of defined length. If the PathLength is chosen to be greater than one, then the individual will overshoot the leader. This path is perturbed randomly.
4.3. Simulated Annealing (SA)
Simulated annealing is one of older algorithm compared to SOMA and DE. It was introduced by Kirkpatrick et al. for the first time [21]. An inspiration for developing this algorithm was annealing of metal. In the process, metal is heated up to temperature near the melting point and then it is cooled very slowly. The purpose is to eliminate unstable particles. In other words, particles are moved towards an optimum energy state. Metal is then in more uniform crystalline structure.
This approach was used in the case of simulated annealing including terms. It starts off from a randomly selected point. Then, a certain number of points (depends on user) are generated in the neighbourhood. The point with the best cost value is selected to be the middle of new neighbourhood (start point for a new loop). However, it is possible to accept also worse value of cost function. The acceptance is based on a probability which decreases with the number of iterations. In the case that the best cost value is in the start point, this one is chosen for the next loop. This approach is basic and some other improvements were done during research in this algorithm.
5. Experimental Results
The main idea is to show that SOMA, DE, and SA are able to solve such problems of symbolic regression–-setting a trajectory–-under analytic programming.
50 simulations were carried out for each algorithm (i.e., 150 simulations in total). SOMA and DE have almost all simulations with positive results; only one case in both algorithms did not reach the extreme. SA was not so successful, only 14 positive results. To show that AP is able to work with arbitrary evolutionary algorithms, we suppose to carry simulations out with genetic algorithms (GAs) and other algorithms, and also parallel computing is intended in this field. Data from all simulations were processed and vizualised in [20, 22].
In simulations made for the purposes of this article, following setting was used to run SOMA, DE, and SA according to Tables 1, 2, and 3, and explanation of each parameter symbol can be found in [15] (SOMA), [20] (DE), and [21] (SA).
Firstly, the results show values of cost function evaluations. This parameter shows good performance of analytic programming. As can be seen in Table 4, the lowest number of cost function evaluations equal 2697 for SA and 3396 for SOMA. DE was also not so far with its 4030 cost function evaluations.
Figure 5 shows the same as Table 4, but in a graphical way, where the diamond means the average value. As can be seen, SA had the lowest average value. However, this might be caused by only 14 cases which were included in the chart while SOMA and DE had 49 positive cases.
Second indicator depicts histogram of successful hits and the number of cost function evaluations for each hit (see Figures 6, 7, and 8). Negative results are not included.
Another creation of histograms can be made from the point of view of number of cases (axe ) which appeared in some interval of cost function values (axe ). This approach can be seen in Figures 9, 10, and 11. Here are all solutions, also bad ones which are represented by higher value than zero.
Next point, which we were interested in, was a number of commands for the ant and number of steps required to eat all baits (Tables 5 and 6). In Table 6, DE found a route which is overcome in the least number of steps can be seen. Sorted lists of pairs, commands and steps, are seen for all 3 algorithms in Table 7. As it is shown, it can be stated that the smallest number of commands does not have to cause the smallest number of steps. Vice versa, the small number of steps does not mean the small set of commands.
Figure 12 depicts that the ant went through all fields; the white “” shows fields which were attended by the ant. The notation (9) contains a set of rules for the ant how to go successfully through the trail. In (10), the whole description of the route can be seen where Ea, So, We, and No mean east, south, west, and north (which cardinal point the ant is turned into). The numbers in brackets are positions on the grid:
6. Conclusions
This contribution deals with an alternative algorithm for symbolic regression. This study shows that this algorithm is suitable not only for mathematical regression but also for setting of optimal trajectory for artificial ant which can be replaced by robots in real world, in industry.
In comparison with standard GP, it can be stated on the basic aforementioned results that AP can solve this kind of problems in shorter times as cost function evaluations are counted.
The aim of this study was not to show that AP is better or worse than GP (or GE when compared), but that AP is also a powerful tool for symbolic regression with support of different evolutionary algorithms.
The main object of this paper was to show that symbolic regression done by AP is able to solve also cases where linguistic terms as, for example, commands for movement of artificial ant or robots in real world are. Here, simulations for 3 algorithms: SOMA, DE, and SA were carried out. As the figures showed, SOMA and DE were more successful in positive results than SA was. This proved that a good performance of AP depends on a choice of suitable robust and powerful evolutionary algorithms.
During simulations carried in this problem following results were reached:
(I)50 simulations for each algorithm means 150 in total for all 3 algorithms.(II)Positive results:(i)49 from 50 simulations for SOMA,(ii)49 from 50 for DE,(iii)and 14 from 50 for SA,which accomplished the required tasks thus analytic programming is able to solve such kind of problems in symbolic regression. This result also says that the basic version of simulated annealing used here is not so powerful tool as other two evolutionary algorithms are. It is supposed that the cost function is very complicated with quite a lot of local optima and, therefore, the simulated annealing was not so successful as SOMA or DE were.(III)Solutions which fulfil conditions which were laid down by Koza [1], concerned to the number of steps, were found (2 by SOMA and 3 by DE). It means 5 solutions were successful under the 400 steps. Moreover, 17 (SOMA) + 20 (DE) + 6 (SA), in total 43 from 150 were successful under the 545 steps which was introduced by Koza [1, 22] as an optimal one.
Future research is key activity in this field. The following steps are to finished simulations with GA and other evolutionary algorithms and to try some other class of problems to show that analytic programming is powerful tool as genetic programming or grammatical evolution are.
Acknowledgments
This work was supported by Grant no. MSM 7088352101 of the Ministry of Education of the Czech Republic and by grants of the Grant Agency of the Czech Republic GACR 102/09/1680.