Ants have always fascinated human beings. What particularly strikes the occasional observer as well as the scientist is the high degree of societal organization that these insects can achieve in spite of very limited individual capabilities (Dorigo et al., 2000). Ants have inspired also a number of optimization algorithms. These algorithms are increasingly successful among researches in computer science and operational research (Blum, 2005; Cordón et al., 2002; Dorigo & Stützle, 2004).
A particular successful metaheuristic—Ant Colony Optimization (ACO)—as a common framework for the existing applications and algorithmic variants of a variety of ant algorithms has been proposed in the early nineties by Marco Dorigo (Dorigo, 1992). ACO takes inspiration from the foraging behavior of some ant species. These ants deposit pheromone on the ground in order to mark some favorable path that should be followed by other members of the colony. ACO exploits a similar mechanism for solving combinatorial optimization problems.
In recent years ACO algorithms have been applied to more challenging and complex problem domains. One such domain is continuous optimization. However, a direct application of the ACO for solving continuous optimization problem is difficult.
The first algorithm designed for continuous function optimization was continuous ant colony optimization (Bilchev & Parmee, 1995) which comprises two levels: global and local; it uses the ant colony framework to perform local searches, whereas global search is handled by a genetic algorithm. Up to now, there are few other adaptations of ACO algorithm to continuous optimization problems: continuous interacting ant colony (Dréo & Siarry, 2002), ACO for continuous and mixed-variable (Socha, 2004), aggregation pheromone system (Tsutsui, 2006), and multilevel ant-stigmergy algorithm (Korošec & Šilc, 2008).
In this chapter we will present so-called Differential Ant-Stigmergy Algorithm (DASA), a new approach to the continuous optimization problem (Korošec, 2006). We start with the DASA description followed by three case studies which show real-world application of the proposed optimization approach. Finally, we conclude with discussion of the obtained results.
2. A differential ant-stigmergy approach
2.1. Continuous optimization
The general continuous optimization problem is to find a set of parameter values, that minimizes a function, of D real variables, i.e.,
To solve this problem, we created a fine-grained discrete form of continuous domain. With it we were able to represent this problem as a graph. This enabled us to use ant-based approach for solving numerical optimization problems.
2.2. The fine-grained discrete form of continuous domainLet be the current value of the i-th parameter. During the searching for the optimal parameter value, the new value, is assigned to the i-th parameter as follows:
Here, is the so-called parameter difference and is chosen from the set
Therefore, for each parameter the parameter difference, has a range from to where b is the so-called discrete base,
With the parameter the maximum precision of the parameter is set. The precision is limited by the computer's floating-point arithmetic. To enable a more flexible movement over the search space, the weight is added to Eq. 2:
2.3. Graph representation
From all the sets where D represents the number of parameters, a so-called search graph with a set of vertices, V, and a set of edges, E, between the vertices is constructed (see Fig. 2). Each set is represented by the set of vertices,
Then we have that
is equal to
Each vertex of the set is connected to all the vertices that belong to the set Therefore, this is a directed graph, where each path from start vertex to any of the ending vertices is of equal length and can be defined with as:
The optimization task is to find a path such that where is currently the best solution, and (using Eq. 9). Additionally, if the objective function is smaller than then the values are replaced with values.
2.4. The differential ant stigmergy algorithm
The optimization consists of an iterative improvement of the currently best solution, by constructing an appropriate path that uses Eq. 9 and returns a new best solution. This is done as follows (see Fig. 3):
Step 1. A solution is manually set or randomly chosen.
Step 2. A search graph is created and an initial amount of pheromone, is deposited on all the vertices from the set according to the Cauchy probability density function
where is the location offset and is the scale factor. For an initial pheromone distribution the standard Cauchy distribution ( and ) is used and each parameter vertices are equidistantly arranged between
Step 3. There are ants in a colony, all of which begin simultaneously from the start vertex. Ants use a probability rule to determine which vertex will be chosen next. More specifically, ant in step moves from a vertex in set to vertex with a probability, prob, given by:
where is the amount of pheromone on vertex
The ants repeat this action until they reach the ending vertex. For each ant, path is constructed. If for some predetermined number of tries we get the search process is reset by randomly choosing new and pheromone re-initialization. New solution is constructed (see Eq. 9) and evaluated with a calculation of
The current best solution, out of solutions is compared to the temporary best solution If is better than then values are replaced with values. In this case is increased (in our case for 1 %) and pheromone amount is redistributed according to the associated path Furthermore, if new is better then then values are replaced with values. So, global best solution is stored. If no better solution is found is decreased (in our case for 3 %).
Step 4. Pheromone dispersion is defined by some predetermined percentage The probability density function is changed in the following way:and
Pheromone dispersion has a similar effect as pheromone evaporation in classical ACO algorithm.
Step 5. The whole procedure is then repeated until some ending condition is met.
It is a well known that ant-based algorithms have problems with convergence. This happens when on each step of the walk there is a large number of possible vertices from which ant can choose from. But this is not the case with the DASA where Cauchy distribution of pheromone over each parameter was used. Namely, such distribution reduces the width of the graph to only few dominant parameter values (i.e., vertices). On the other hand, with proper selection of the discrete base, b, we can also improve the algorithm's convergence (larger b reduces the search graph size).
2.5. Software implementation and parameter settings
The DASA was implemented in the Borland® Delphi™ programming language (see Fig. 4). The computer platform used to perform the optimizations was based on AMD Opteron™ 2.6-GHz processors, 2 GB of RAM, and the Microsoft® Windows® XP 32-bit operating system.
The DASA parameters were set as follows: the number of ants, was set to 10, the pheromone evaporation factor, was set to 0.2, and the maximum parameter precision, dependent on the discrete step of each parameter.
3. Case studies
Case studies presented in this section are related to the development of a new dry vacuum cleaner (Korošec et al., 2007; Tušar et al., 2007). In particular, the efficiency of the turbo-compressor unit (see Fig. 5) was improved.
3.1. An electric motor power losses minimization
Home appliances, such as vacuum cleaners and mixers, are generally powered by a universal electric motor (UM). These appliances need as low as energy consumption, that is, input power, as possible, while still satisfying the needs of the user by providing sufficient output power. The optimization task is to find the geometrical parameter values that will generate the rotor and the stator geometries resulting in the minimum power losses.
There are several invariable and variable parameters that define the rotor and stator geometries. The invariable parameters are fixed; they cannot be altered, either for technical reasons (e.g., the air gap) or because of the physical constraints on the motor (e.g., the radius of the rotor's shaft). The variable parameters do not have predefined optimum values. Some variable parameters are mutually independent and without any constraints. Others are dependent, either on some invariable parameters or on mutually independent ones.
In our case, 11 mutually independent variable parameters defining the rotor and stator geometries are subject to optimization (see Fig. 6): rotor yoke thickness, rotor external radius, rotor pole width, stator yoke horizontal thickness, stator yoke vertical thickness, stator middle-part length, stator internal edge radius, stator teeth radius, stator slot radius, stator width, and stator height, We optimized only 10 parameters (the ratio between and was set constant) of the UM rotor and stator geometries.
Power losses, are the main factor affecting the efficiency of a UM. The efficiency of a UM, is defined as the ratio of the output power, to the input power,
and it is very dependent on various power losses (Sen, 1997):
The overall copper losses, occurring in the rotor and the stator slots are as follows:
where stands for each slot, is the current density, is the slot area, is the copper's specific resistance and is the length of the winding turn. The calculation of the iron losses, is less exact because iron has non-linear magnetic characteristics. The iron losses are therefore separated into two components: the hysteresis losses and the eddy-current losses. This means the motor's iron losses can be expressed with the formula
where is the eddy-current material constant at 50Hz, is the hysteresis material constant at 50Hz, is the maximum magnetic flux density, is the frequency of the magnetic field density in the rotor, is the frequency of the magnetic field density in the stator, is the mass of the rotor, and is the mass of the stator. Three additional types of losses occurring in the UM, i.e., the brush losses, the ventilation losses, and the friction losses, depend mainly on the speed of the motor. When optimizing the geometries of the rotor and the stator, the motor's speed is assumed to be fixed; hence and have no impact on the motor's efficiency, and so these losses are not significantly affected by the geometries of the rotor and the stator. Therefore, in our case, the overall copper and iron losses can be used to estimate the cost function:
Our goal is to find such parameter-value settings, where is minimized.
The DASA starts its search with the existent solution. The stopping criterion was set to 1,400 calculations. The calculation of a single solution via the ANSYS Multiphysics simulation takes approximately two minutes, which means that the execution of 1,400 calculations needs about two days. The optimization method was run 20 times. The obtained results in terms of the UM power losses are presented statistically in Table 1.
The engineering rotor and stator design (the existent solution) results in power losses of 177.9 W, and can be seen in Fig. 7 (left). The figure shows the magnetic flux density in the laminations (higher magnetic flux density, causes higher power losses). Figure 7 (right) present a typical example of a feasible rotor and stator geometry, with power losses of 129.1 W. This solution has very low iron losses in the rotor due to its small size and in spite of its high magnetic saturation. The small rotor and its saturation are compensated by large stator poles that ensure large enough a magnetic flux. This design is completely feasible from the technical and production points of view.
3.2. An electric motor casing stiffness maximization
The casing is a part of the dry vacuum cleaner motor which is built into vacuum cleaners. The casing is basically an axisymmetric shell structure which is built of steel suitable for forming. For this procedure which consists of eleven different phases it is important that the radii are growing or do not change while the height is growing. That is an important rule which must not be broken during the geometry optimization. The goal of optimization is to preserve the stiffness while using a thinner shell structure and consequently save material and reduce costs.
The computational model of the casing comprises 26 different parameters which generate the geometry variations. The classical parameters are: radius on the top, radius on the side, radius at the groove for brushes, deviation of the hole for diffuser fixation, radius of roundness at the beginning of the bearing groove, radius on the top of the bearing groove, radius on the top of air culverts, radius at the bottom of air culverts, roundness of air culverts, angle of slope, angle of culverts span, slope of bearing groove, height of bearing groove, and slope of the brushes groove,
Besides the above listed 14 classical parameters there are 12 more (p 10 to p 21), with which the ribs on areas A, B and C are defined (Fig. 8). These are essential to improve the stiffness so we will search for those which maximize the cost function. With these parameters we can generate a lot of combinations of different rib shapes, depending on the number of given paths and accuracy of intervals for points deviations, positions and number of ribs. In our case there were several billion combinations for ribs. It is obvious that instead of a defined path at position j, we could define displacements for individual points but this would result in many more parameters and time consuming calculation procedure. The question still remains what would be the benefit of that. From previous designs there exist some reasonable shapes for ribs which have been defined intuitively and these then change their position, magnitude and number, which are probably sufficient for a good result, as there has been over ten different rib shapes defined earlier.
Besides the parameters for the creation of ribs there are other classical parameters which also have an influence on their shape. Namely, if the radii or are varying then this is reflected in the shape of ribs of model A or if the height of bearing groove, and slope of groove, are modified then this will be reflected in the ribs of model C, etc.
The cost function is relatively simply defined when loadings are deterministic. Usually we optimize such cases with the minimization of displacements or stresses, maximization of elasticity, etc. (Dražumerič & Kosel, 2005). It is difficult to determine the cost function for a dynamically loaded assembly where the loads are stochastic. So the question is how to define the stiffness of a shell structure if we do not know the lateral loads. Let us consider a certain arbitrarily shaped plane shell. We can write its potential energy, as:
If we consider that if the Poissons shear modulus is neglected, the stress tensor, and strain tensor, are the following:
Here is displacement in direction and is elasticity modulus. The potential energy can then be written as:
We can see that in this case the potential energy is an integral of squares of stresses over the volume. We define also the kinetic energy, for plate vibrations:
where represents the density of material. Here the second derivative of displacement is approximated with respect to time with the product of squared displacement and eigen frequency,
For conservative systems the maximal potential energy is equal to maximal kinetic energy and from this follows the expression for cost function which can be in our case the stiffness itself:
The DASA was run 20 times and each run consisted of 2,000 calculations. The cost function was calculated by the ANSYS Multiphysics simulation tool. The obtained results are presented statistically in Table 2.
The results of optimization were quite surprising (Fig. 9 bottom row). It was expected that the ribs would form on all the given surface but they did not. The ribs were more distinct where the vertical part of the surface was turning into the horizontal one (ribs A). There were no ribs at the surface where the stator is in a tight fit, probably because of pre-stressing of the casing through the stator. Also in the groove for brushes (ribs B) there were no distinct ribs, probably because the groove itself is a kind of rib. Instead of this there is a given slope for a groove which was not there before. The radii of roundings were in most cases bigger than before as we can clearly see at the air culverts (Fig. 9 middle bottom).
3.3. A turbo-compressor aerodynamic power maximization
Radial air impellers are the basic components of many turbo-machines. In the following we will concentrate on relatively small impellers and subsonic speeds used in a dry vacuum cleaner. Our main aim was to find an impeller shape that has a higher efficiency, i.e., greater aerodynamic power, than the one currently used in production.
An impeller is constructed from blades, an upper and a lower side. The sides enclose the blades and keep them together. The blades, which are all the same, were the main part of the optimization. The geometry of a blade is shown in Fig. 10, where the gray color represents the blade. The method of modeling is as follows: we construct the points at specific locations, draw the splines through them and spread the area on the splines. Once a blade is made an air channel must be constructed in a similar way.
In Fig. 10a the point 1 has two parameters: the radius and the angle Similarly, the points 2, 5 and 6 have parameter pairs and and and The points 3 and 4 are fixed on the axis. This is because the impeller must have a constant outer radius and the outer side of the blade must be parallel to the axis. On the other hand, the outer angle of the blade and the angle of the spline at points 3 and 4, can be varied. Analogously, the angles and are the inner-blade angles for the upper and lower edges of the blade at the input, respectively.
In Fig. 10c the points 1, 2, and 3 form the upper spline, and the points 4, 5, and 6, the lower spline. Between the points 1 and 6 is the point 7, which defines the spline of the input shape of the blade. In this figure, the points 1, 2, 5, and 6 have the parameters and respectively, describing their heights.
Point 3 stays on the axis and point 4 has a constant height In other words, the designer of the impeller must know at least the outer diameter and the height The parameters and describe the input angles of the lower and upper parts of the blade with respect to the – plane. Similarly, the parameters and describe the outer blade angle with respect to the same plane.
In Fig. 10b the meaning of point 7 is explained more precisely. The parameters and define the radius, height, and angle, respectively. The radius and angle dictate where the point should appear with respect to the – plane and the height with respect to the – plane. Similarly, the angles and are needed to define the starting and ending angles of the spline constructed between the points 1, 7, and 6.
If we look closely at Fig. 10c then we can see the contour surrounding the blade. This is the air channel with the following parameters: the inner radius (see Fig. 10a), which is needed for the hexahedral mesh (explained later), the air intake radius the air outflow radius the bolt radius the bolt height and the impeller height
In this way we have successfully modeled the impeller geometry with 34 parameters. For each parameter we have a predefined search interval with a given discrete step. Therefore, the size of the search space can be obtained as the product of the number of possible settings over all the parameters. It turns out that there are approximately 3 1034 possible solutions.
The mesh of the air channel between the two blades is constructed with more than 6,000 hexahedral elements. The boundary conditions are zero velocity at all the solid regions and symmetry boundary conditions at the fluid regions. At the influx and outflux the intake velocity, and reference pressure, are defined, respectively. The intake velocity is parabolically distributed, because we expect that the intake flow is laminar and so:
Here, is a velocity dependent on the stream, which further depends on time, as we shall see later, is the upper radius, defined before, and is the radius within the limits from to The reference pressure, is set to zero.
With respect to the maximum time, , the flux is:
where is the influx area and is the current time.
The distribution of the relative pressure can be used to estimate the cost function. The average pressure, is chosen from the air-intake area. Finally, the aerodynamic power, which represents the cost function, is as follows:
where is the output pressure at the radius and l/s is the flux near the desired optimum performance of the impeller. Our goal is to find such parameter-value settings, where is maximized.
The DASA was run 10 times and each run consisted of 2,000 CFD calculations. For the CFD calculations we used the ANSYS Multiphysics package. A single CFD calculation takes approximately seven minutes. The obtained results, in terms of aerodynamic power, are presented statistically in Table 3.
Results show that we were able to increase aerodynamic power for approximately 20 %. Figure 11 shows a 3D view of the existent and two optimized impellers (best and worst of 10 runs).
In this chapter the Differential Ant-Stigmergy Algorithm (DASA) was presented, where the main goal was an evaluation of the DASA on some real-world engineering applications. Case studies were selected from a R&D project where more efficient turbo-compressor for dry vacuum cleaner was developed. Here the DASA was used to improve the efficiency of an electric motor, increase casing stiffness, and increase impeller's aerodynamic power. In all these cases the improvement was evident.
In (Korošec & Šilc, 2009a) we also compared the DASA to the eleven state-of-the-art stochastic algorithms for continuous optimization which were presented in the CEC'2005 Special Session on Real Parameter Optimization. The algorithms are:
BLX-GL50, a hybrid real-coded genetic algorithm with female and male differentiation (García-Martínez & Lozano, 2005),
BLX-MA, a real-coded memetic algorithm with adaptive local search parameters (Molina et al., 2005),
CoEVO, an evolutionary algorithm with mutation step co-evolution (Pošík, 2005),
DE, a differential evolution (Rönkkönen et al., 2005),
EDA, a continuous estimation of distribution algorithm (Yuan & Gallagher, 2005),
FEA, a flexible evolutionary algorithm (Alonso et al., 2005),
G-CMA-ES, a restart covariance matrix adaptation evolutionary strategy with increasing population size ( Auger & Hansen, 2005a ),
G-DE, a guided differential evolution (Bui et al., 2005),
K-PCX, a population-based steady-state procedure (Sinha et al., 2005),
L-CMA-ES, an advanced local search evolutionary algorithm (Auger & Hansen, 2005b), and
SPC-PNX, a steady-state real-parameter genetic algorithm (Ballester et al., 2005).
Obtained results confirmed that the DASA is comparable to above algorithms and therefore generally applicable to global optimization problems.
The experimental results have shown that the DASA ranks among the best algorithms for real-life-like applications. Its main advantage is in consistent problem solving which is indicated by 19 rankings in top third, 4 ranking in middle third and only 2 in bottom third.
Statistical analysis following the Bonferroni-Dunn post-hoc test with showed that the DASA is significantly better than 8 out of 11 compared algorithms.
The algorithm was also applied to dynamic optimization problems with continuous variables proposed for CEC’2009 Special Session on Evolutionary Computation in Dynamic and Uncertain Environments (Korošec & Šilc, 2009b). If we compare the DASA to:
the results show that the DASA can find reasonable solutions for all of the problems. It can be seen that the DASA performs not many worse than a self-adaptive differential evolution and much better than the other three algorithms. One obvious advantage is that was no need any changes to the original algorithm. So, it can be used as such for both cases of numerical optimization, static and dynamic. Furthermore, the algorithm is unsusceptible to different types of changes and can be used with very limited knowledge about problem, only maximal dimension and input problem parameters.