Versatility of Simulated Annealing with Crystallization Heuristic: Its Application to a Great Assortment of Problems

This chapter is related to several aspects of optimization problems in engineering. Engineers usually mathematically model a problem and create a function that must be minimized, like cost, required time, wasted material, etc. Eventually, the function must be maximized. This function has different names in the literature: objective function, cost function, etc. We will refer to it in the chapter as objective function. There is a wide range of possibilities for the problems and they can be classified in different ways. At first, the values of the parameters can be continuous, discrete (integers), cyclic (angles), intervals, and combinatorial. The result of the objective function can be continuous, discrete (integers) or intervals. One very difficult class of problems have continuous parameters and discrete objective function, this type of objective function has very weak sensibility. This chapter shows the versatility of the simulated annealing showing that it can have different possibilities of parameters and objective functions.


Introduction
Optimization problems are common in engineering, and the problem must be optimized regarding one function (usually called as cost function). The cost function has several parameters used to model the optimization problem. Several methods have been proposed in the literature to optimize functions, and the method is selected according the cost function and parameters characteristics.
Deterministic methods usually have an initial point, called as seed. Using the gradient of the cost function, the optimization method starts on the initial point and converges to the optimal parameter value. The parameters are usually continuous [1].

Simulated annealing with crystallization heuristic
The SA is an optimization algorithm that uses a random local search and is able to find the global optimum solution. It was proposed by [5,14], and it was inspired by the metal annealing process, which starts with a high temperature and reduces the temperature in small increments until reaching the minimum temperature. In each step, the new solution is accepted if it improves the current solution. Otherwise, it is only accepted if a probability factor P T ð Þ is more than a random number [10,14,15]. The P T ð Þ can be determined by where T is the current temperature, and ΔE is the the energy state variation which is the difference between the objective function F x ð Þ of the new candidate a and the current solution b, as presented in ΔE ¼ F a ð Þ À F b ð Þ. Algorithm 1 describes a possible implementation of a conventional SA. It initializes with a random solution and an initial temperature, which is an important parameter. There are two loops. The external loop controls the temperature. The temperature decreases according to a geometric cooling schedule with factor α, and it stops when the frozen state is reached. The cooling schedule, in this case the value of α, must be carefully chosen, as the convergence to the optimal distribution strongly depends on it. In the majority of the applications explained here, we used α ¼ 0:98. The internal loop, performs the random local search, until the thermal equilibrium is reached. The random search consists of modifying a single parameter and a new candidate x * is obtained. The new objective function value is evaluated and compared with the cost of the current solution. Then, it is either accepted or rejected according to (1). The convergence conditions are usually controlled by algorithm parameters. They influence the quality of the final result as well as the speed of the algorithm.
< initial temperature > 4: while < Global condition not satisfied > do 5: while < Local condition not satisfied > do 6: x * < modify a parameter > 7: The random search in the SA algorithm can use different strategies to improve finding a new solution. One of them is the SA with crystallization heuristic proposed by [9,16]. The SA has two phases when the search is performed: exploration and refinement. Figure 1 shows the two phases and their connection with the crystallization factor. The exploration of the domain space, usually happens in higher temperatures. The refinement of the solution happens in lower temperatures. There is no clear interface between the two phases.
For problems with combinatorial and integer parameters, the generation of the new candidate can be easily performed. In the case of combinatorial parameters, the algorithm exchanges the position of two parameters. In the case of integer parameters, the algorithm increments, decrements or keeps the current value. There are possibilities, considering that integer parameters have lower and upper bounds, it is possible to generate a new value performing a random selection in the range.
The control of continuous parameters is challenging. The crystallization heuristic proposed by [9,16] considers a fixed maximum step, and the SA controls the probability density. Figure 1 shows that wider probability density allows larger jumps with higher probability. While, thinner probability density allows smaller jumps with higher probability. However, a larger jump still possible to happen with very small probability. This feature allows the SA to escape from local minima.
Each continuous parameter has a crystallization factor, which is represented by c j . Considering the crystallization heuristic, the new candidate is determined by summing c j times a random number, this procedure determines a Bates distribution. It is represented by where x is the current solution, x * is the new candidate, they are represented by vectors. Δr j is the fixed step size associated with continuous parameter j. e j represents the selected parameter, e j is a vector with all elements equal to zero except one, the position associated with the j-th parameter. A common value assigned to Δr j depends on the search interval, as Δr j ¼ max j À min j À Á =4. This value will provide enough exploration for the algorithm.
It remains to explain the procedure in which the SA controls the crystallization factor for each parameter. The SA modifies only one parameter at a time, and according the decision of accepting or rejecting the new candidate, an action is performed. The procedure is represented in Figure 2. If the new candidate is rejected, it is assumed that the SA is performing exploration and to increase accepted solutions the crystallization factor associated with this parameter is increased. The increase in the crystallization factor will reduce the probability of larger jumps for this parameter. On the other hand, if the new candidate is accepted, it is assumed that the SA is performing refinement and to increase exploration the crystallization factor associated with this parameter is reduced. The reduction of the crystallization factor will enlarge the probability of larger jumps for this parameter.  x < random initial solution > 3: T 0 < initial temperature > 4: Δr < range of the parameters > 5: while < Global condition not satisfied > do 6: while < Local condition not satisfied > do 7: j < select parameter to modify > 8: x if ΔE < 0 then 11: x x * 12: if c j > 1 then 13: T i T iÀ1 * α 28: end while Algorithm 2 describes the implementation of SA with crystallization heuristic. In this algorithm, the crystallization heuristic adjusts the modification to increase the acceptance of new solutions. The crystallization factor c is initialized with 1 for all continuous parameters, and the fixed step size Δr i is defined as 25% of the search range. For example, if the search happens between À100 and þ100, in this case Δr i ¼ 50.

Application in cutting and packing
The field of operations research concentrates many real world applications of optimization techniques in engineering, as it essentially e aims to increase the efficiency of industry operations [8,9,[16][17][18][19]. Interest in this area has grown recently due to advances in optimization algorithms and the importance of the waste and pollution reduction, which is usually an effect of an improved solution.
Among operations research subjects, cutting and packing (C&P) problems can be highlighted due to its importance and singular combination of geometry and optimization. These characterises are even more prominent in irregular packing problems, which involves simple polygonal items. Essentially, C&P problems consists of assigning a set of small items to a set of set of large containers, maintaining some geometric restrictions, while minimizing an objective function.
In this section, a SA based solution for the irregular packing problem is proposed. It adopts a discrete objective function, as it is related to the number of items in the container, while using continuous parameter for items rotations.

Irregular bin packing problem with free rotations
In the 2D irregular single bin packing problem, given a collection of items, one must place a subset of these inside a rectangular container with the aim of minimizing the unused space inside the bin. Each item is represented by a simple polygon and may be rotated by any angle. The main geometric restrictions dictates that no two items may overlap and there should be no protrusion from the container.
Consider a set of items P ¼ P 1 , P 2 , … , P n f gand a rectangular container C. The layout can be represented by the translation vector t 1 , t 2 , … , t n f g , a rotation vector r 1 , r 2 , … , r i f gand an assignment set T , which contains the subset items placed inside the container. The irregular bin packing problem can be described as the minimization of the difference between the area of the container A C ð Þ and the assigned set of items A T ð Þ, as where P r ð Þ represents an item rotate by r, operator i P ð Þ expresses only the interior of P and P ⊕ t indicates a translation t applied to P.

Solution using SA
One of the main challenges in irregular packing is the complexity of managing the geometric constraint, which results in fewer proposed solutions in the literature [20]. In order to optimize the layout without compromising the geometric feasibility, two main strategies are often employed. The first is to define a constructive heuristic, which places one item at a time, usually at the bottom-left position. In this approach, the optimization algorithm only controls the placement order; a popular solution is to employ genetic algorithms [2,21]. The alternative is to allow the items to move freely inside the container while applying penalization on the item overlaps [22][23][24].
In [9], however, a different approach was adopted: a SA based algorithm was employed to directly control the placement order, as well as the position and rotation of each item. Items were placed sequentially, using the parameters given by the SA, until no more items fitted the container. Then, the objective function, which was the unoccupied area of the container, which is directly related to the subset of placed items, was evaluated.
The main difficulty was the definition of the item position parameter, which should always correspond to a valid placement, i.e., without overlap or protrusion from the container. It was given by the collision free region (CFR), a polygon describing the allowed placement region, as shown in Figure 3. The CFR can be obtained using modified Boolean operations on polygons, described in [25]. A continuous parameter evaluated to a position along the perimeter of the CFR, and the closest vertex was chosen as the placement position of the item. The rotation was defined by a second controlled variable, and it had to be applied prior to the determination of the CFR.
Therefore, the solution optimization basically consisted of a SA controlling the placement and rotation parameter of each item in the layout. At each iteration, one parameter for a single item was changed, or two items were swapped in the placement order. Then, the cost was evaluated and the new solution was accepted according to (1). The crystallization factor described in Section 2 was applied to the rotation parameters.

Results and discussion
Six broken glass puzzle instances were created to evaluate the performance of the algorithm. These instances have a known optimal solution, which are displayed in Figure 4. Therefore, the effectiveness of the algorithm can be measure by its success rate, which relates to the number of executions converging to the optimal solution.
The tests were executed on a Phenom 9550 2.21GHz and the convergence condition for the SA was that (1) the cost variation for the final temperature was zero, and (2) the final solution has the lowest cost found. The initial temperature was adjusted by admitting initial 50% acceptance and a geometric cooling schedule was adopted with α ¼ 0:98 (as shown in Algorithm 1). Table 1 shows the results data obtained for 30 executions of each instance. Results show that only very simple problems were solved optimally in every executions. This indicates that the SA algorithm should be complemented with heuristics approaches to increase the convergence rate. Nevertheless, given the complexity of the problem with free rotations, the fact that it found the optimal solutions for all instances is important, as some could not be achieved by enforcing simple heuristics such as bottom-left or larger first.
One important characteristic of the irregular bin packing problem with free rotations is that the objective function is discrete and some of the parameters are continuous. This is illustrated in Figure 5, which shows discrete cost values for different values for the leftmost item rotation. The SA solution was not affected by

Application in TO
TO is a mathematical approach to determine distribution of material in a design domain such that the performance becomes maximized. The definition of performance could be different in each application. The physic of problem and desired application determines the objective function and the constraints. Depending on the problem, different methods have been developed in the literature [26,27].
For the TO problems with well defined objective functions and constraints, the gradient based algorithms have been widely used [28][29][30]. The sensitivity information converges the final results to the optimized topology. This method shows fast convergence and low computational costs. The gradient based TO methods use Optimality Criteria (OC), Method of Moving Asymptotes (MMA), Sequential Linear Programming (SLP), etc. to optimize the objective functions [31,32]. In the cases that the objective function or its derivatives are not mathematically modeled or hard to calculate, using the non-gradient based algorithms are more advantages [33,34]. In the non-gradient based TO methods, GA and SA are the two most popular algorithms of optimization [3,[35][36][37]. Even these methods have high computational costs, but can optimize the topology without need to calculation of the derivatives and sensitivity information.
Using the SA method in the non-gradient based TO is beneficial because of reaching global minimum as well as provide the information of convergence.  The information of convergence such as the number of accepted and rejected solutions could be used in the evaluation of TO results. The available SA for TO uses random search to generate new solutions, while by using SA with crystallization heuristic [10], the search for new accepted solutions has been improved. After finishing the optimization process, a density filter can reduce gray area and discontinued regions.

SA for non-gradient based TO
The structural TO to minimize compliance in beams is a classic problem. This problem has been solved with gradient based methods in the literature [38] and the results have been used to verify TO with crystallization heuristic SA (as described in Section 2). The problem of minimizing compliance can be represented as the problem of minimizing strain energy, modeled as where F , U, and K are, respectively, external force, elastic deformation, and stiffness. The constraint is volume fraction that represents the final optimized topology should have less or equal volume of the design domain. By discretizing the design domain to N square elements, the total strain energy can be calculated as where x e is the density of each element varying from a minimum value (to avoid singularity in matrix calculation) to 1. p is the penallization parameter, u e is elastic deformation for element e and k e is the stiffness for element e. The penalization factor p penalizes the intermediate gray area in the Solid Isotropic Material with Penalization (SIMP) method to reduce gray area. In this case, the TO has continuous parameters and it is solved using SA with crystallization heuristic, as described in Section 2.
A new heuristic is included, after reaching the thermal equilibrium, the domain is regularized by filtering. The new density of each element after filtering gets some effects from the adjacent elements, as where w e is the weighting function which is the filter radius minus distance of each adjacent element. It should be noted that the weighting function is zero outside of the filter radius and the density changes just inside the filer domain. The design domain and loading conditions for cantilever and half-MBB beam problems are shown in Figure 6. A comparison of obtained compliance from proposed method and some results from the literature is shown in Table 2.
As shown in Table 2, the results from this non-gradient based TO method are very close to the gradient based results. The main advantage of the proposed method is that there is no need to calculate derivatives of the objective function.

SA for multi-objective TO
The TO objective function can be complex and combine two or more objective functions. In such situations, the solution is not necessarily unique and comes with a set of optimum solutions called Pareto Front. The Pareto Front curve shows the solutions where none of the objective functions can be improved without degrading the other objective functions. The Pareto Front curve can be used for trade-off the suitable solution within this set instead of considering the full range of every parameter. The traditional TO usually optimizes one objective function while considering the other objective functions fixed or as a constraint [39]. The SA showed also this ability to incorporate multiple objective functions, like the one presented by the CoAnnealing [40] and AMOSA [41].   The CoAnnealing was used to solve TO problems and the Pareto Front was obtained; particularly, compliance and volume fraction as considered as cost functions in a cantilever beam, as shown in Figure 7. The points on the curve showed a good agreement with the results for that volume fraction and compliance [42]. The results shows that Coannealing can be used for the multiple objective TO problems.

Application to curve interpolation
The curve fitting is a CAD (Computer Aided Design) study area. This problem consists in determining a curve from a set of points, as shown in Figure 8. Two different types of curves can be determined, an interpolating curve that passes through the set of points and an approximating curve that passes near to the points. The set of points is shown in blue points, the interpolating curve in red line and the approximating curve the dashed line.
In this section, a SA approach is proposed to determine an approximation curve from a sequence of points, in which the control points are the continuous parameters and index corresponding points of the sequence are discrete parameters. All these parameters are adjusted by the SA. This study was started by Ueda et al. [43]. The piece-wise Bézier curve is used to overcome some difficulties, as explained in the following section.

Piece-wise Bézier curve
There are several curves that can be used in the curve fitting problem, each one of them has features that is a drawback. The Bézier curve's control points modify all the curve, and regions that the curve is already fitted. On the other hand, the control points of the B-spline have a local influence, i.e., just a region of the curve is modified. One feature of the Bézier curve that aids in the fitting problem is that the curve always interpolates the first and last control point, while in the B-spline these points usually are not interpolated. It is possible to interpolates these points in the B-spline. However, a higher number of parameter is needed to achieve such feature, compared to the Bézier curve.
A piece-wise Bézier curve overcome the problem of the global influence of the control point. This curve is a sequence of cubic Bézier curves as shown in Figure 9, in which the last control point of curve p 3 is the first control point of the following curve. The determination of the second control point of the second curve p 4 is given by with β being a proportional factor that ensures the weak-G1 continuity between the curves, i.e., the tangent vector of the end of the first curve has the same direction but not necessary the same intensity of the tangent of the start of the second curve. Ueda et al. [44] proposed an algorithm for automatic evaluate the number of piece-wise Bézier curves necessary to interpolate a sequence of points.

SA for the approximation curve determination
Using a piece-wise Bézier curve has a advantage in two aspects. First, the curve is divided in several pieces, in which each curve segment approximates only part of the sequence. Second, it is known that a Bézier curve always interpolates its end points. The used cost function is with p i are the control points defining the curve, d is the sequence of points, P v d k ð Þ is the projection of point d k in the curve P, L P u ð Þ ð Þis the length of the curve, P is a piece-wise Bézier curve and W is a weight factor. This object function consists of two parts. P mÀ1 k¼1 d k À P v d k ð Þ k k 2 is the first part and it represents the discrepancy between the piece-wise curve and the sequence of points. The first and last points in the sequence are not considered as they already are on the determined piece-wise Bézier curves. For every point in the sequence d k it is determined its projection in the piece-wise Bézier curve P v d k ð Þ [45]. L P u ð Þ ð Þis the second part, representing a regularization. It is used to avoid the over-fitting problem, the algorithm will search for short and smooth curves (and avoid long and sinuous curves).
SA needs to adjust the parameters defining a piece-wise Bézier curve that minimizes (8). Each piece-wise Bézier curve has two internal control points (as shown in Figure 9). The coordinates for each of these control points are continuous parameters and they are controlled using the crystallization heuristic. The connection point between two piece-wise Bézier curves (as point p 3 shown in Figure 9) is represented as an integer parameter. It is an index representing a point from the given sequence.
Consider the example from Figure 10, the continuous parameters are the coordinate of control points p 1 , p 2 and p 5 . Control point p 4 is determined by modifying β from (7). Control points p 0 and p 6 are the first and last point of the sequence of points; and control point p 3 is a connecting point. The definition of this point is the discrete parameter adjusted by the SA.
The continuous parameters are controlled by the crystallization heuristic. The discrete parameters have specific operator, a random value is picked from a uniform distribution between 0, 1 ½ , if this value is greater than 0:5 the discrete parameter is increased by 1, otherwise it is decreased by 1.

Results and discussion
A curve is defined and sampled, and a random noise is added to each sampled point creating an artificial example to be tested by the proposed algorithm. Two example were tested, shown in Figure 11. Figure 11(a) is sampled with 101 points and Figure 11(b) is sampled with 300 points. The initial temperature is set as 100 and for each temperature there are 1000 iteration or 500 accepted solutions, and an adaptive cooling schedule [45] was adopted with a variable α. A minimum of 5 accepted solutions is used as stop criteria. Table 3 show the result of these two test examples. The objective function is higher in the second example, once N pts is higher as well the curve is longer. N iter is close in both example, there are less than 7% of difference between both tests.
The proposed interpolation algorithm was applied in different applications showing its versatility and robustness. Ueda et al. [46] proposed an algorithm to determine the curves which interpolates the boundary of a point cloud in space.  Ueda et al. [47] proposed an algorithm to determine the defect zone from a 3D scanned turbine blade. Ueda et al. [48] also proposed an algorithm to determine open boundaries in point clouds with symmetry (like injuries in a skull).

Conclusions
This chapter describes the SA with crystallization heuristic and three different applications: cutting and packing, topology optimization and curve interpolation. The cutting and packing problem determines the layout for a set of items to be placed inside a container with fixed dimensions. The cost function is the unused area inside the container. The items can be translated and rotated in the process of determining the layout. This cost function is discrete and the parameters are continuous. As there is no need of any gradient information, it can optimize discrete cost functions.
The topology optimization problem has continuous parameters. The number of parameters is much larger when compared with the cutting and packing problem. Usually, the cost function is one of the two possibilities: volume fraction and compliance. As there are two possibilities for cost function, it is considered the application of the CoAnnealing, which is a multi-objective version of the SA with crystallization heuristic. The last application is the curve interpolation. It has two different types of parameters: continuous and integer. The objective function is composed of two parts: the error between the the given points and the interpolating curve, and the length of the curve. The length of the curve is a regularization to force the curve to be shorter and smoother. The curve interpolation algorithm was used in different applications.
These examples show that the SA with crystallization heuristic is very generic, versatile and easy to implement.