Open access peer-reviewed chapter - ONLINE FIRST

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

Written By

Tiago G. Goto, Hossein R. Najafabadi, Guilherme C. Duran, Edson K. Ueda, André K. Sato, Thiago C. Martins, Rogério Y. Takimoto, Hossein Gohari, Ahmad Barari and Marcos S.G. Tsuzuki

Reviewed: May 25th, 2021 Published: July 6th, 2021

DOI: 10.5772/intechopen.98562

Engineering Problems - Uncertainties, Constraints and Optimization Techniques Edited by Marcos Sales Guerra Tsuzuki

From the Edited Volume

Engineering Problems - Uncertainties, Constraints and Optimization Techniques [Working Title]

Dr. Marcos Sales Guerra Tsuzuki and Prof. Rehab O. Abdel Rahman

Chapter metrics overview

22 Chapter Downloads

View Full Metrics


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.


  • Simulated Annealing
  • Optimization Problems
  • Cutting Packing
  • Electrical Impedance Tomography
  • Topological Optimization
  • Curve Interpolation

1. 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]. However, engineering problems have problems which cannot be solved by this approach.

Some cost functions have several minimum, and the convergence strongly depends on the seed (or initial point). To overcome this issue, several meta-heuristics were proposed, like Genetic Algorithms (GA) [2, 3], Particle Swarm Optimization (PSO) [4], Simulated Annealing (SA) [5], and others. Optimization methods based on meta-heuristics do not need the gradient information. For some problems, it is very hard to determine the gradient. This is particularly true for problems with discrete cost function. The existing meta-heuristics, usually, are focused on some specific type of parameters. For example, PSO can be applied to engineering problems with continuous parameters. GA can be applied to engineering problems with integer or fixed precision parameters.

SA was first proposed to solve combinatorial problems [5]. Later on, Corana et al. [6] extended the SA to incorporate continuous parameters. It was shown that the proposal made by Corana et al. [6] only converged to the global optimum in specific problems. Ingber [7] improved the SA to overcome this issue, however he used the cost function gradient for some calculations. Martins and Tsuzuki [8] proposed the SA with crystallization heuristic which does not use any gradient information.

This chapter shows the versatility of the SA with crystallization heuristic. Initially, the SA with crystallization heuristic is explained in Section 2. It was used to solve different problems related to Cutting and Packing, particularly the case with discrete cost function and cyclic continuous parameters [8, 9, 10] will be explained in Section 3. A third application is Topological Optimization (TO), it has a large number of parameters, like the EIT problem. This application is explained in Section 4. The curve interpolation problem, explained in Section 5, has two different types of parameters: continuous and integer. It is relevant to point that the SA with crystallization heuristic has been applied with success to interval based objective function evaluation [11, 12, 13]. The conclusions are in Section 6.


2. 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 PTis more than a random number [10, 14, 15]. The PTcan be determined by


where Tis the current temperature, and ΔEis the the energy state variation which is the difference between the objective function Fxof the new candidate aand the current solution b, as presented in ΔE=FaFb.

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 xis 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.

Algorithm 1: Conventional SA

 1: i0

 2: x<random initial solution >

 3: T0<initial temperature >

 4: while<Global condition not satisfied>do

 5:   while<Local condition not satisfied>do

 6:    x<modify a parameter >

 7:    ΔE=FxFx

 8:    ifΔE<0or random01<eΔE/Tithen

 9:      xx

10:    end if

11:   end while

12:   ii+1

13:   TiTi1α

14: end while

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.

Figure 1.

The crystallization heuristic represented by SA controlling each parameter crystallization factor.

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 cj. Considering the crystallization heuristic, the new candidate is determined by summing cjtimes a random number, this procedure determines a Bates distribution. It is represented by


where xis the current solution, xis the new candidate, they are represented by vectors. Δrjis the fixed step size associated with continuous parameter j. ejrepresents the selected parameter, ejis a vector with all elements equal to zero except one, the position associated with the j-th parameter. A common value assigned to Δrjdepends on the search interval, as Δrj=maxjminj/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.

Figure 2.

During the optimization process, the SA has phases: exploration and refinement. In the exploration phase, the parameters have low crystallization and can perform larger jumps. The inverse happens in the refinement phase where the parameters have high crystallization.

Algorithm 2: SA with Crystallization Heuristic.

 1: i0; c111

 2: x<random initial solution >

 3: T0<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:    xx+1cj1cjrandom1212Δrjej

 9:    ΔE=FxFx

10:    ifΔE<0then

11:     xx

12:     ifcj>1then

13:      cjcj1            ▷ Positive Feedback

14:     end if

15:    else

16:     ifrandom01<eΔE/Tithen

17:       xx

18:       ifcj>1then

19:         cjcj1          ▷ Positive Feedback

20:       end if

21:     else

22:       cjcj+1           ▷ Negative Feedback

23:     end if

24:    end if

25:  end while

26:  ii+1

27:  TiTi1α

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 cis initialized with 1for all continuous parameters, and the fixed step size Δriis defined as 25%of the search range. For example, if the search happens between 100and +100, in this case Δri=50.


3. 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.

3.1 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=P1P2Pnand a rectangular container C. The layout can be represented by the translation vector t1t2tn, a rotation vector r1r2riand 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 ACand the assigned set of items AT, as


where Prrepresents an item rotate by r, operator iPexpresses only the interior of Pand Ptindicates a translation tapplied to P.

3.2 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.

Figure 3.

Example of collision free region for itemP.

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.

3.3 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.

Figure 4.

Instances for the irregular bin packing problem.

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.

LF Fails54,42610.6163.6%
BL Fails51170.16100.0%

Table 1.

The results for the SA based irregular bin packing solution. Nitems: total number of items. Nconv: average number of iterations. Tconv: average time untill convergence (in seconds). Pconv: success rate.

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 this issue, as the results show, it can handle both continuous and discrete parameters and objective function. Moreover, the continuous rotation convergence for each item was improved by adopting the crystallization factor.

Figure 5.

Discrete cost behaviour with single continuous parameter variation.


4. 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.

4.1 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 Kare, 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 Nsquare elements, the total strain energy can be calculated as


where xeis the density of each element varying from a minimum value (to avoid singularity in matrix calculation) to 1. pis the penallization parameter, ueis elastic deformation for element eand keis the stiffness for element e. The penalization factor ppenalizes 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 weis 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.

Figure 6.

The design domain and loading for cantilever beam (left) and half-MBB beam (right).

Volume fractionCantilever in [38]Cantilever in SAHalf-MBB in [38]Half-MBB in SA

Table 2.

The results for compliance of cantilever and half-MBB beam for different volume fractions.

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.

4.2 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.

Figure 7.

The Pareto Front curve for minimization of compliance and weight in cantilever problem.


5. Application to curve fitting

The problem of curve fitting is an essential Computer Aided Design (CAD) problem with applications in various engineering fields including but not limited to digital metrology, robotics, path planning, and data modeling. The problem consists of constructing a curve from a set of discrete points, as shown in Figure 8. The two main types of curve fitting can be defined, as an approximating curve fitting, when the constructed curve approximates the location of the points in the dataset, and an interpolating curve fitting, when the curve passes exactly through the set of points. The curve approximation has application in many engineering problems when a certain level of uncertainty exists in the dataset. This always can be expected in data points collected by an experimental process, when sources of uncertainties including the equipment errors, environmental effects, human errors, measurement resolution, etc., represent a combined level of uncertainty [43]. Due to the non-systematic nature of the uncertainties in the data, the approximating curve fitting process aims to recognise the true pattern in the data points instead of exactly passing through them [44]. The three major computation tasks to complete the approximating curve fitting include Point Measurement Planning (PMP), Substitute Geometry Estimation (SGE), and Deviation Zone Evaluation (DZE). Reducing or controlling the level of uncertainty in the constructed curve have been studied comprehensively at PMP by proper selection of the datapoints [45]. It is also addressed in SGE by improving the curve or surface fitting algorithms, using an enhanced optimization processes to avoid trapping in local minima, and by using iterative fitting approaches with monitoring some indicators for the level of uncertainty [46]. Various approaches have been also presented to measure and monitor the level of uncertainty typically at DZE stage by modeling the pattern and nature of the resulting deviations of the processed data points from the approximated curves or surfaces [47].

Figure 8.

Approximating curve (dashed) versus the interpolating curve (red).

On the contrary, the curve interpolation is applicable when the datapoints are known to be fairly accurate and are conducted with no accommodation for any level of uncertainty in the datapoints. Figure 8 presents a set of datapoints with its corresponding approximating and interpolating curves. The set of datapoints is shown by blue dots, the interpolating curve is presented by solid red line, and the approximating curve is presented by black dashed line.

Since a certain level of uncertainty typically exists in the most of the engineering problems, developing methodologies to control or reduce the level of uncertainty in the finally constructed curve is highly important. In this section, a SA approach is proposed to determine an approximation curve from a sequence of points. In this approach, the control points are continuous parameters and index corresponding points of the sequence are discrete parameters. All these parameters are adjusted by SA. This study was started by Ueda et al. [48]. The developed methodology employs piece-wise Bézier curve structure to solve the problem, as explained in the followings.

5.1 Piece-wise Bézier curve

There are several curve structures that can be used for the purpose of solving a curve fitting problem. However, each one has it own advantages and drawbacks. In a Bézier curve, the control points influence the entire curve globally including the regions that the curve is already fitted. On the other hand, the control points in a B-spline have only a local influence, i.e., by changing a group of control points only a certain region of the curve is modified. One feature of the Bézier curve structure that is beneficial in the presented fitting approach here is that the resulting Bézier curve always interpolates the first and last control point, while in the B-spline curve structure, these points usually are not interpolated. It is possible to interpolates these points in a B-spline. However, a higher number of optimization parameters 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 p3is the first control point of the following curve. The determination of the second control point of the second curve p4is given by

Figure 9.

Piece-wise cubic Bézier curve with2curve segments, pointsp0,p1,p2andp3defines the first curve and pointsp3,p4,p5andp6defines the second one.


with βbeing a proportional factor that ensures the weak-G1continuity 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. [49] proposed an algorithm for automatic evaluate the number of piece-wise Bézier curves necessary to interpolate a sequence of points.

5.2 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 piare the control points defining the curve, dis the sequence of points, Pvdkis the projection of point dkin the curve P, LPuis the length of the curve, Pis a piece-wise Bézier curve and Wis a weight factor. This object function consists of two parts. k=1m1dkPvdk2is 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 dkit is determined its projection in the piece-wise Bézier curve Pvdk[50]. LPuis 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 p3shown 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 p1, p2and p5. Control point p4is determined by modifying βfrom (7). Control points p0and p6are the first and last point of the sequence of points; and control point p3is a connecting point. The definition of this point is the discrete parameter adjusted by the SA.

Figure 10.

Each Bézier curve segment approximates just one part of the sequence of points.

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 01, if this value is greater than 0.5the discrete parameter is increased by 1, otherwise it is decreased by 1.

5.3 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 101points and Figure 11(b) is sampled with 300points. The initial temperature is set as 100and for each temperature there are 1000iteration or 500accepted solutions, and an adaptive cooling schedule [50] was adopted with a variable α. A minimum of 5accepted solutions is used as stop criteria.

Figure 11.

Two artificial created examples tested by the proposed algorithm. Black points are the noisy sampled points, each curve segment is in a different color with its respective control points represented as squares.

Table 3 show the result of these two test examples. The objective function is higher in the second example, once Nptsis higher as well the curve is longer. Niteris close in both example, there are less than 7%of difference between both tests.


Table 3.

The results for the SA based curve interpolation solution. Npts: total number of points. NSegm: number of curve segments. Niter: number of iterations. Cost: objective function.

The proposed interpolation algorithm was applied in different applications showing its versatility and robustness. Ueda et al. [51] proposed an algorithm to determine the curves which interpolates the boundary of a point cloud in space. Ueda et al. [52] proposed an algorithm to determine the defect zone from a 3D scanned turbine blade. Ueda et al. [53] also proposed an algorithm to determine open boundaries in point clouds with symmetry (like injuries in a skull).


6. 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.



G. C. Duran is supported by CNPq (process 140.299/2020-3). H. R. Najafabadi is supported by FAPESP (grant 2019/03453-2). M. S. G Tsuzuki is partially supported by CNPq (process 311.195/2019-9). A. K. Sato is supported by FUSP/Petrobras.




Temperature Geometric Cooling Parameter


Function Area


Proportional Factor to Ensure Weak-G1 continuity


Rectangular Conteiner


Crystallization Factor


Energy Variation


Step Variation for Variable j


Position Vector for Variable j(all elements are zero, except one)


External Force


Objective Function


Stiffness for Element e


Interior of the Item


Length of Curve Pu




Number of Square Elements


Penalization Parameter


Control Point


Set of Items




Item Rotated by r


Item Translated by t


Piece-wise Bézier Curve


Projection of Point dkin Curve Pu

Space of Continuous Numbers




Strain Energy




Assigned Set of Items


Temperature at iteration i


Elastic Deformation


Elastic Deformation for Element e


Current Solution


Solution Candidate


Density of Element e


Weighting Function


Weight Used in the Approximation Curve




Archived Multiobjective Simulated Annealing


Computer-Aided Design


Collision Free Region


Cutting & Packing


Electrical Impedance Tomography


Finite Element Method


Genetic Algorithm


Maximum a Posteriori


Method of Moving Asymptotes


Optimality Criteria


Probability Density Function


Particle Swarm Optimization


Random Sample Consensus


Simulated Annealing


Solid Isotropic Material with Penalization


Sequential Linear Programming


Topology Optimization


  1. 1. J. V. Burke and M. C. Ferris, “A Gauss-Newton method for convex composite optimization,”Math Program, vol. 71, no. 2, pp. 179–194, 1995.
  2. 2. P. R. Pinheiro, B. Amaro Júnior, and R. D. Saraiva, “A random-key genetic algorithm for solving the nesting problem,”Int J Comp Integ M, vol. 29, no. 11, pp. 1159–1165, 2016.
  3. 3. S. Wang, K. Tai, and M. Y. Wang, “An enhanced genetic algorithm for structural topology optimization,”Int J Numer Meth Eng, vol. 65, no. 1, pp. 18–44, 2006.
  4. 4. E. G. Castro and M. S. G. Tsuzuki, “Swarm intelligence applied in synthesis of hunting strategies in a three-dimensional environment,”Expert Syst Appl, vol. 34, no. 3, pp. 1995–2003, 2008.
  5. 5. S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,”Science, vol. 220, no. 4598, pp. 671–680, 1983.
  6. 6. A. Corana, M. Marchesi,C. martini, and S. Ridella, “Minimizing multimodal functions of continuous variables with the simulated annealing algorithm,”ACM Transactions on Mathematical Software, vol. 13, pp. 262–280, 1987.
  7. 7. L. Ingber, “Very fast simulated re-annealing,”Mathematical and Computer Modelling, vol. 12, no. 8, pp. 967–973, 1989.
  8. 8. T. C. Martins and M. S. G. Tsuzuki, “Simulated annealing applied to the rotational polygon packing,”IFAC Proc Vol (IFAC-PapersOnline), vol. 12, no. PART 1, 2006.
  9. 9. T. C. Martins and M. S. G. Tsuzuki, “Rotational placement of irregular polygons over containers with fixed dimensions using simulated annealing and no-fit polygons,”J Braz Soc Mech Sci, vol. 30, no. 3, pp. 205–212, 2008.
  10. 10. T. C. Martins, A. K. Sato, and M. S. G. Tsuzuki, “Adaptive neighborhood heuristics for simulated annealing over continuous variables,” inSimulated Annealing - Advances, Applications and Hybridizations(M. S. G. Tsuzuki, ed.), pp. 3–20, INTECHOpen, 2012.
  11. 11. T. C. Martins, E. D. L. B. Camargo, R. G. Lima, M. B. P. Amato, and M. S. G. Tsuzuki, “Electrical impedance tomography reconstruction through simulated annealing with incomplete evaluation of the objective function,” in33rd IEEE EMBC, pp. 7033–7036, 2011.
  12. 12. T. C. Martins and M. S. G. Tsuzuki, “Electrical impedance tomography reconstruction through simulated annealing with total least square error as objective function,” in34th IEEE EMBC, pp. 1518–1521, 2012.
  13. 13. R. S. Tavares, T. C. Martins, and M. S. G. Tsuzuki, “Electrical impedance tomography reconstruction through simulated annealing using a new outside-in heuristic and GPU parallelization,”J Phys: Conf Series, vol. 407, p. 012015, dec 2012.
  14. 14. V. Černỳ, “Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm,”J Optimiz Theory App, vol. 45, no. 1, pp. 41–51, 1985.
  15. 15. E. H. L. Aarts and P. J. M. van Laarhoven, “Statistical cooling: a general approach to combinatorial optimization problems,”Philips J Res, vol. 40, no. 4, pp. 193–226, 1985.
  16. 16. T. C. Martins and M. S. G. Tsuzuki, “Placement over containers with fixed dimensions solved with adaptive neighborhood simulated annealing,”B Pol Acad Sci-Tech, vol. 57, pp. 273–280, 2009.
  17. 17. A. K. Sato, T. C. Martins, and M. S. G. Tsuzuki, “Rotational placement using simulated annealing and collision free region,”IFAC Proceedings Volumes, vol. 43, no. 4, pp. 234–239, 2010.
  18. 18. A. K. Sato, T. C. Martins, and M. S. G. Tsuzuki, “An algorithm for the strip packing problem using collision free region and exact fitting placement,”CAD, vol. 44, no. 8, pp. 766–777, 2012.
  19. 19. A. K. Sato, T. de Castro Martins, and M. S. G. Tsuzuki, “A pairwise exact placement algorithm for the irregular nesting problem,”Int J Comp Integ M, vol. 29, no. 11, pp. 1177–1189, 2016.
  20. 20. J. A. Bennell and J. F. Oliveira, “The geometry of nesting problems: A tutorial,”Eur J Oper Res, vol. 184, no. 2, pp. 397 – 415, 2008.
  21. 21. L. R. Mundim, M. Andretta, and T. A. de Queiroz, “A biased random key genetic algorithm for open dimension nesting problems using no-fit raster,”Expert Syst Appl, vol. 81, pp. 358–371, 2017.
  22. 22. J. Egeblad, B. K. Nielsen, and A. Odgaard, “Fast neighborhood search for two- and three–dimensional nesting problems,”Eur J Oper Res, vol. 183, pp. 1249–1266, 2007.
  23. 23. A. Elkeran, “A new approach for sheet nesting problem using guided cuckoo search and pairwise clustering,”Eur J Oper Res, vol. 231, no. 3, pp. 757–769, 2013.
  24. 24. A. K. Sato, T. C. Martins, A. M. Gomes, and M. S. G. Tsuzuki, “Raster penetration map applied to the irregular packing problem,”Eur J Oper Res, vol. 279, no. 2, pp. 657–671, 2019.
  25. 25. A. K. Sato, T. C. Martins, and M. S. G. Tsuzuki, “Collision free region determination by modified polygonal boolean operations,”CAD, vol. 45, no. 7, pp. 1029–1041, 2013.
  26. 26. M. P. Bendsoe and O. Sigmund,Topology optimization: theory, methods, and applications. Springer Science & Business Media, 2013.
  27. 27. T. R. Marchesi, R. D. Lahuerta, E. C. N. Silva, M. S. G. Tsuzuki, T. C. Martins, A. Barari, and I. Wood, “Topologically optimized diesel engine support manufactured with additive manufacturing,”IFAC-PapersOnLine, vol. 28, no. 3, pp. 2333–2338, 2015.
  28. 28. O. Sigmund and K. Maute, “Topology optimization approaches,”Struct Multidiscip O, vol. 48, no. 6, pp. 1031–1055, 2013.
  29. 29. G. I. Rozvany, “A critical review of established methods of structural topology optimization,”Struct Multidiscip O, vol. 37, no. 3, pp. 217–237, 2009.
  30. 30. A. Shukla and A. Misra, “Review of optimality criterion approach scope, limitation and development in topology optimization,”Int J Adv Eng Techn, vol. 6, no. 4, p. 1886, 2013.
  31. 31. S. Rojas-Labanda and M. Stolpe, “Benchmarking optimization solvers for structural topology optimization,”Struct Multidiscip O, vol. 52, no. 3, pp. 527–547, 2015.
  32. 32. F. A. Gomes and T. A. Senne, “An slp algorithm and its application to topology optimization,”Comput Appl Math, vol. 30, no. 1, pp. 53–89, 2011.
  33. 33. O. Sigmund, “On the usefulness of non-gradient approaches in topology optimization,”Struct Multidiscip O, vol. 43, no. 5, pp. 589–596, 2011.
  34. 34. D. Guirguis, W. W. Melek, and M. F. Aly, “High-resolution non-gradient topology optimization,”J Comput Phys, vol. 372, pp. 107–125, 2018.
  35. 35. N. Garcia-Lopez, M. Sanchez-Silva, A. Medaglia, and A. Chateauneuf, “A hybrid topology optimization methodology combining simulated annealing and simp,”Comput Struct, vol. 89, no. 15-16, pp. 1512–1522, 2011.
  36. 36. G. Y. Cui, K. Tai, and B. P. Wang, “Topology optimization for maximum natural frequency using simulated annealing and morphological representation,”AIAA J, vol. 40, no. 3, pp. 586–589, 2002.
  37. 37. M. Ohsaki, “Genetic algorithm for topology optimization of trusses,”Comput Struct, vol. 57, no. 2, pp. 219–225, 1995.
  38. 38. O. Sigmund, “A 99 line topology optimization code written in matlab,”Struct Multidiscip O, vol. 21, no. 2, pp. 120–127, 2001.
  39. 39. K. Suresh, “A 199-line matlab code for pareto-optimal tracing in topology optimization,”Struct Multidiscip O, vol. 42, no. 5, pp. 665–679, 2010.
  40. 40. T. C. Martins and M. S. G. Tsuzuki, “EIT image regularization by a new multi-objective simulated annealing algorithm,” inIEEE EMBC 2015, pp. 4069–4072, IEEE, 2015.
  41. 41. T. C. Martins, A. V. Fernandes, and M. S. G. Tsuzuki, “Image reconstruction by electrical impedance tomography using multi-objective simulated annealing,” inIEEE ISBI 2014, pp. 185–188, jul 2014.
  42. 42. H. R. Najafabadi, T. G. Goto, T. C. Martins, A. Barari, and M. S. G. Tsuzuki, “Multi-objective topology optimization using simulated annealing method,” inICGG 2020, pp. 343–353, Springer.
  43. 43. A. Barari, “Automotive body inspection uncertainty associated with computational processes,”Int J Veh Des, vol. 57. no. 2-3, pp. 230-241, 2011.
  44. 44. A. Barari, H. A. ElMaraghy, and P. Orban, “NURBS representation of estimated surfaces resulting from machining errors,”Int J Comp Integ M, vol. 22, no. 5, pp. 395-410, 2009.
  45. 45. A. Lalehpour, C. Berry, and A. Barari, “Adaptive data reduction with neighbourhood search approach in coordinate measurement of planar surfaces,”J M Syst, vol. 45. pp. 28-47, 2017.
  46. 46. A. Barari, and S. Mordo, “Effect of sampling strategy on uncertainty and precision of flatness inspection studied by dynamic minimum deviation zone evaluation,”J Metrol Qual Eng, vol. 4, no. 1, pp. 3-8, 2013.
  47. 47. A. Lalehpour, and A. Barari, “Developing skin model in coordinate metrology using a finite element method,”Measurement, vol. 109, pp. 149-159, 2017.
  48. 48. E. K. Ueda, M. S. G. Tsuzuki, R. Y. Takimoto, A. K. Sato, T. C. Martins, P. E. Miyagi, and R. S. U. Rosso Jr, “Piecewise Bézier curve fitting by multiobjective simulated annealing,”IFAC-PapersOnLine, vol. 49, pp. 49–54, jan 2016.
  49. 49. E. K. Ueda, T. C. Martins, and M. S. G. Tsuzuki, “Planar curve fitting by simulated annealing with feature points determination,”IFAC-PapersOnLine, vol. 51, pp. 290–295, jan 2018.
  50. 50. E. K. Ueda, A. K. Sato, T. C. Martins, R. Y. Takimoto, R. S. U. Rosso Jr, and M. S. G. Tsuzuki, “Curve approximation by adaptive neighborhood simulated annealing and piecewise Bézier curves,”Soft Comput, vol. 24, p. 18821–18839, 2020.
  51. 51. E. K. Ueda, M. S. Tsuzuki, and A. Barari, “Piecewise Bézier curve fitting of a point cloud boundary by simulated annealing,” inINDUSCON 2018, pp. 1335–1340, jan 2019.
  52. 52. E. K. Ueda, A. Barari, A. K. Sato, and M. S. G. Tsuzuki, “Detection of defected zone using 3D scanning data to repair worn turbine blades,” in21st IFAC World Congress, pp. 10666–10670, 2020.
  53. 53. E. K. Ueda, A. Barari, and M. S. G. Tsuzuki, “Determination of open boundaries in point clouds with symmetry,” inICGG 2020, pp. 332–342, Springer.

Written By

Tiago G. Goto, Hossein R. Najafabadi, Guilherme C. Duran, Edson K. Ueda, André K. Sato, Thiago C. Martins, Rogério Y. Takimoto, Hossein Gohari, Ahmad Barari and Marcos S.G. Tsuzuki

Reviewed: May 25th, 2021 Published: July 6th, 2021