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

Abstract

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.

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 $P\left(T\right)$ is more than a random number [10, 14, 15]. The $P\left(T\right)$ can be determined by

$$P\left(T\right)={exp}^{-\frac{\mathrm{\Delta }E}{T}}$$E1

where $T$ is the current temperature, and $\mathrm{\Delta }E$ is the the energy state variation which is the difference between the objective function $F\left(x\right)$ of the new candidate $a$ and the current solution $b$, as presented in $\mathrm{\Delta }E=F\left(a\right)-F\left(b\right)$.

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 $\alpha$, and it stops when the frozen state is reached. The cooling schedule, in this case the value of $\alpha$, 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 $\alpha =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^{\ast }$ 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.

Algorithm 1: Conventional SA

 1: $i\leftarrow 0$

 2: $x\leftarrow <$ random initial solution $>$

 3: $T_0\leftarrow <$ initial temperature $>$

 4: while$<$ Global condition not satisfied$>$do

 5:   while$<$ Local condition not satisfied$>$do

 6:    $x^{\ast }\leftarrow <$ modify a parameter $>$

 7:    $\mathrm{\Delta }E=F\left({\mathit{x}}^{\ast }\right)-F\left(\mathit{x}\right)$

 8:    if$\mathrm{\Delta }E<0$ or $\mathit{\text{random}}\left(0,1\right)<e^{-\mathrm{\Delta }E/T_i}$then

 9:      $\mathit{x}\leftarrow {\mathit{x}}^{\ast }$

10:    end if

11:   end while

12:   $i\leftarrow i+1$

13:   $T_i\leftarrow T_{i-1}\ast \alpha$

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 $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

$$x^{\ast }=x+\frac{1}{c_j}\sum _{k=1}^{c_j}\text{random}\left(-1/2,1/2\right)\cdot \mathrm{\Delta }{r}_j\cdot e_j.$$E2

where $x$ is the current solution, $x^{\ast }$ is the new candidate, they are represented by vectors. $\mathrm{\Delta }{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 $\mathrm{\Delta }{r}_j$ depends on the search interval, as $\mathrm{\Delta }{r}_j=\left({\mathit{max}}_j-{min}_j\right)/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: $i\leftarrow 0$; $c\leftarrow \left(1,1,\dots ,1\right)$

 2: $x\leftarrow <$ random initial solution $>$

 3: $T_0\leftarrow <$ initial temperature $>$

 4: $\mathrm{\Delta }r\leftarrow <$ range of the parameters $>$

 5: while$<$ Global condition not satisfied$>$do

 6:  while$<$ Local condition not satisfied$>$do

 7:    $j\leftarrow <$ select parameter to modify $>$

 8:    ${\mathit{x}}^{\ast }\leftarrow \mathit{x}+\frac{1}{c_j}{\sum }_1^{c_j}\mathit{\text{random}}\left(-\frac{1}{2},\frac{1}{2}\right)\cdot \mathrm{\Delta }{r}_j\cdot {\mathit{e}}_j$

 9:    $\mathrm{\Delta }E=F\left({\mathit{x}}^{\ast }\right)-F\left(\mathit{x}\right)$

10:    if$\mathrm{\Delta }E<0$then

11:     $\mathit{x}\leftarrow {\mathit{x}}^{\ast }$

12:     if$c_j>1$then

13:      $c_j\leftarrow c_j-1$             ▷ Positive Feedback

14:     end if

15:    else

16:     if$\mathit{\text{random}}\left(0,1\right)<e^{-\mathrm{\Delta }E/T_i}$then

17:       $\mathit{x}\leftarrow {\mathit{x}}^{\ast }$

18:       if$c_j>1$then

19:         $c_j\leftarrow c_j-1$           ▷ Positive Feedback

20:       end if

21:     else

22:       $c_j\leftarrow c_j+1$            ▷ Negative Feedback

23:     end if

24:    end if

25:  end while

26:  $i\leftarrow i+1$

27:  $T_i\leftarrow T_{i-1}\ast \alpha$

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 $\mathrm{\Delta }{r}_i$ is defined as $25\%$ of the search range. For example, if the search happens between $-100$ and $+100$, in this case $\mathrm{\Delta }{r}_i=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 $\mathcal{P}=\left\{P_1,P_2,\dots ,P_n\right\}$ and a rectangular container $\mathcal{C}$. The layout can be represented by the translation vector $\left\{t_1,t_2,\dots ,t_n\right\}$, a rotation vector $\left\{r_1,r_2,\dots ,r_i\right\}$ and an assignment set $\mathcal{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\left(\mathcal{C}\right)$ and the assigned set of items $A\left(\mathcal{T}\right)$, as

$$\begin{array}{ll}minimize& A\left(\mathcal{C}\right)-A\left(\mathcal{T}\right)\\ \text{subject}\text{to}& i\left(P_i\left(r_i\right)\oplus t_i\right)\cap \left(P_j\left(r_j\right)\oplus t_j\right)=\varnothing i,j\in \mathcal{T}\text{and}i<j\\ & \left(P_i\left(r_i\right)\oplus t_i\right)\subseteq \mathcal{C}1\le i\in \mathcal{T}\\ & 0\le r_i<2\pi 1\le i\in \mathcal{T}\\ & t_i\in {\Re }^{2}1\le i\in \mathcal{T}\end{array}$$E3

where $P\left(r\right)$ represents an item rotate by $r$, operator $i\left(P\right)$ expresses only the interior of $P$ and $P\oplus t$ indicates a translation $t$ applied 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 item $P$.

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 $\alpha =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.

Instance	Nitems	Nconv	$T_{}$conv(s)	$P_{}$conv
Tangram	7	130,169	722.83	58.8%
Hole	7	12,847	28.89	93.5%
Simple	4	24,604	49.58	100.0%
Concave	5	17,912	43.96	86.7%
LF Fails	5	4,426	10.61	63.6%
BL Fails	5	117	0.16	100.0%

Table 1.

The results for the SA based irregular bin packing solution. $N_{\mathit{\text{items}}}$: total number of items. $N_{\mathit{\text{conv}}}$: average number of iterations. $T_{\mathit{\text{conv}}}$: average time untill convergence (in seconds). $P_{\mathit{\text{conv}}}$: 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

$$minimizeS=U^T\mathit{KU},\text{subjected}\text{to}:\mathit{KU}=\mathcal{F}$$E4

where $\mathcal{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

$$S=\sum _{e=1}^N{\left(x_e\right)}^p{u}_e^T{k}_e{u}_e$$E5

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

$$x_{\text{filter}}=\frac{{\sum }_{e=1}^N{w}_e{x}_e}{{\sum }_{e=1}^N{w}_e}$$E6

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.

Figure 6.

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

Volume fraction	Cantilever in [38]	Cantilever in SA	Half-MBB in [38]	Half-MBB in SA
0.5	75.05	79.87	82.37	89.65
0.6	62.89	65.70	69.19	73.36
0.7	55.27	57.10	60.61	63.39
0.8	49.89	51.29	54.85	56.63
0.9	46.29	47.36	50.57	51.79

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 ${\mathbf{p}}_3$ is the first control point of the following curve. The determination of the second control point of the second curve ${\mathbf{p}}_4$ is given by

Figure 9.

Piece-wise cubic Bézier curve with $2$ curve segments, points ${\mathit{p}}_0,{\mathit{p}}_1,{\mathit{p}}_2$ and ${\mathit{p}}_3$ defines the first curve and points ${\mathit{p}}_3,{\mathit{p}}_4,{\mathit{p}}_5$ and ${\mathit{p}}_6$ defines the second one.

$${\mathbf{p}}_4={\mathbf{p}}_3-\beta \cdot \left({\mathbf{p}}_2-{\mathbf{p}}_3\right),$$E7

with $\beta$ 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. [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

$$F\left(\left\{{\mathbf{p}}_i\right\}\right)=W\cdot \sum _{k=1}^{m-1}{‖{\mathbf{d}}_k-{\mathbf{P}}_v\left({\mathbf{d}}_k\right)‖}^2+\left(1-W\right)\cdot L\left(\mathbf{P}\left(u\right)\right).$$E8

with ${\mathbf{p}}_i$ are the control points defining the curve, $\mathbf{d}$ is the sequence of points, ${\mathbf{P}}_v\left({\mathbf{d}}_k\right)$ is the projection of point ${\mathbf{d}}_k$ in the curve $\mathbf{P}$, $L\left(\mathbf{P}\left(u\right)\right)$ is the length of the curve, $\mathbf{P}$ is a piece-wise Bézier curve and $W$ is a weight factor. This object function consists of two parts. ${\sum }_{k=1}^{m-1}{‖{\mathbf{d}}_k-{\mathbf{P}}_v\left({\mathbf{d}}_k\right)‖}^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 ${\mathbf{d}}_k$ it is determined its projection in the piece-wise Bézier curve ${\mathbf{P}}_v\left({\mathbf{d}}_k\right)$ [50]. $L\left(\mathbf{P}\left(u\right)\right)$ 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 ${\mathbf{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 ${\mathbf{p}}_1$, ${\mathbf{p}}_2$ and ${\mathbf{p}}_5$. Control point ${\mathbf{p}}_4$ is determined by modifying $\beta$ from (7). Control points ${\mathbf{p}}_0$ and ${\mathbf{p}}_6$ are the first and last point of the sequence of points; and control point ${\mathbf{p}}_3$ is 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 $\left[0,1\right]$, if this value is greater than $0.5$ the 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 $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 [50] was adopted with a variable $\alpha$. A minimum of $5$ accepted 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 $N_{\mathit{pts}}$ is higher as well the curve is longer. $N_{\mathit{\text{iter}}}$ is close in both example, there are less than $7\%$ of difference between both tests.

Example	Npts	Nsegm	Niter	Cost
1	101	4	508,404	10.37
2	300	5	473,680	134.97

Table 3.

The results for the SA based curve interpolation solution. $N_{\mathit{pts}}$: total number of points. $N_{\mathit{\text{Segm}}}$: number of curve segments. $N_{\mathit{\text{iter}}}$: number of iterations. $\mathit{\text{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.

Acknowledgments

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.

Nomenclature

$\alpha$

Temperature Geometric Cooling Parameter

$A\left(\right)$

Function Area

$\beta$

Proportional Factor to Ensure Weak-G1 continuity

$\mathcal{C}$

Rectangular Conteiner

$C_j$

Crystallization Factor

$\mathrm{\Delta }E$

Energy Variation

$\mathrm{\Delta }{r}_j$

Step Variation for Variable $j$

$e_j$

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

$F_e$

External Force

$F\left(x\right)$

Objective Function

$k_e$

Stiffness for Element $e$

$i\left(P_k\right)$

Interior of the Item

$L\left(\mathbf{P}\left(u\right)\right)$

Length of Curve $\mathbf{P}\left(u\right)$

$K$

Stiffness

$N$

Number of Square Elements

$p$

Penalization Parameter

${\mathbf{p}}_k$

Control Point

$\mathcal{P}$

Set of Items

$P_k$

Item

$P_k\left(r\right)$

Item Rotated by $r$

$P_k\oplus t$

Item Translated by $t$

$\mathbf{P}\left(u\right)$

Piece-wise Bézier Curve

${\mathbf{P}}_v\left({\mathbf{d}}_k\right)$

Projection of Point ${\mathbf{d}}_k$ in Curve $\mathbf{P}\left(u\right)$

$\Re$

Space of Continuous Numbers

$r$

Rotation

$S$

Strain Energy

$t$

Translation

$\mathcal{T}$

Assigned Set of Items

$T_i$

Temperature at iteration $i$

$U$

Elastic Deformation

$u_e$

Elastic Deformation for Element $e$

$x$

Current Solution

$x^{\ast }$

Solution Candidate

$x_e$

Density of Element $e$

$w_e$

Weighting Function

$W$

Weight Used in the Approximation Curve

Abbreviations

AMOSA

Archived Multiobjective Simulated Annealing

CAD

Computer-Aided Design

CFR

Collision Free Region

C&P

Cutting & Packing

EIT

Electrical Impedance Tomography

FEM

Finite Element Method

GA

Genetic Algorithm

MAP

Maximum a Posteriori

MMA

Method of Moving Asymptotes

OC

Optimality Criteria

PDF

Probability Density Function

PSO

Particle Swarm Optimization

RANSAC

Random Sample Consensus

SA

Simulated Annealing

SIMP

Solid Isotropic Material with Penalization

SLP

Sequential Linear Programming

TO

Topology Optimization