Assembly sequence planning of rigid and flexible parts

Journal of Manufacturing Systems 36 (2015) 128–146

Contents lists available at ScienceDirect

Journal of Manufacturing Systems journal homepage: www.elsevier.com/locate/jmansys

Technical Paper

Assembly sequence planning of rigid and flexible parts Somayé Ghandi, Ellips Masehian ∗ Faculty of Engineering, Tarbiat Modares University, Tehran 14115-143, Iran

a r t i c l e

i n f o

Article history: Received 20 November 2014 Received in revised form 29 April 2015 Accepted 19 May 2015 Keywords: Assembly sequence planning (ASP) Flexible assembly parts Assembly stress matrix (ASM) Scatter search Parameter tuning TOPSIS method

a b s t r a c t Assembly sequence planning (ASP) is the process of computing a sequence of assembly motions for constituent parts of an assembled final product. ASP is proven to be NP-hard and thus its effective and efficient solution has been a challenge for the researchers in the field. Despite the fact that most assembled products like ships, aircrafts and automobiles are composed of rigid and flexible parts, no work exists for assembly/disassembly sequence planning of flexible parts. This paper lays out a theoretical ground for modeling the deformability of flexible assembly parts by introducing the concept of Assembly stress matrix (ASM) to describe interference relations between parts of an assembly and the amount of compressive stress needed for assembling flexible parts. Also, the Scatter Search (SS) optimization algorithm is customized for this problem to produce high-quality solutions by simultaneously minimizing both the maximum applied stress exerted for performing assembly operations and the number of assembly direction changes. The parameters of this algorithm are tuned by a TOPSIS-Taguchi based tuning method. A number of ASP problems with rigid and flexible parts were solved by the presented SS and other algorithms like Genetic and Memetic algorithms, Simulated Annealing, Breakout Local Search, Iterated Local Search, and Multistart Local Search, and the results and their in-depth statistical analyses showed that the SS outperformed other algorithms by producing the best-known or optimal solutions with highest success rates. © 2015 The Society of Manufacturing Engineers. Published by Elsevier Ltd. All rights reserved.

1. Introduction The Assembly Sequence Planning (ASP) problem concerns with finding a sequence of collision-free operations o1 , . . ., on that bring the assembly parts p1 , . . ., pn together, having given the geometry of the final product A and the positions of parts in the final product. ASP plays a key role in the whole lifecycle of a product and has a great impact on variation propagation, production quality, and efficiency of the assembly process [29]. Assembly of manufactured goods constitutes over half of the total production time in manufacturing industry, and about one-third to half of labor costs [32]. The complexity of finding an optimal sequence increases linearly with the size of the space of all potential assembly sequences in the case of exhaustive search. In fact, the ASP problem is classified as an NP-hard problem and thus cannot be solved in polynomial time for large instances as the total number of potential sequences is given by the permutations of parts, n! [22]. In recent years, intensive research efforts have been put in developing intelligent methods to solve the ASP problem, as they have been able to improve the efficiency of finding optimal assembly

∗ Corresponding author. Tel.: +98 2182884939. E-mail address: [email protected] (E. Masehian).

sequences while avoiding the combinatorial explosion problem as the number of assembly components increases. After Bonneville et al. [4] implemented Genetic Algorithm (GA) in ASP for the first time in 1995, many soft computing/metaheuristic algorithms were developed to solve the problem, such as Simulated Annealing (SA) [25,17], Artificial Immune System (AIS) [5], GA, Ant Colony Optimization (ACO) [37,15], Artificial Neural Networks (ANN) [6,16,31,7], Memetic Algorithm (MA) [11], enhanced Harmony Search (HS) [38], Imperialist Competitive Algorithm (ICA) [43], Particle Swarm Optimization (PSO) [42,36,22,39], Firefly Algorithm (FA) [41], Psychoclonal algorithm [34], and Bacterial Chemotaxis algorithm [44]. An inclusive survey on the ASP and its methods is presented in [18], which overviews the elements of sequence planning such as finding a feasible sequence, determining an optimal sequence according to one or more operational criteria, representing the space of feasible assembly sequences in different ways, applying search and optimization algorithms, and satisfying precedence constraints existing between subassemblies. Also, a survey of works that have applied particularly soft computing approaches in ASP has been presented in [27], covering the years 2001 to 2011. In all of the abovementioned works, a main assumption is that all parts of the assembly are rigid (not deformable) and their shapes do not change during the assembly process. This assumption is a

http://dx.doi.org/10.1016/j.jmsy.2015.05.002 0278-6125/© 2015 The Society of Manufacturing Engineers. Published by Elsevier Ltd. All rights reserved.

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

simplifying one since flexible parts introduce additional degrees of complexity to the ASP problem. In fact, most real-world assembled products like ships, aircrafts, and automobiles are composed of rigid and flexible parts, and so automatic generation of assembly sequence plans for such products requires the deformability of flexible parts to be taken into account. The solution to the ASP problem of a product with rigid and flexible parts depends heavily on the model that is used to simulate the deformation behavior of flexible parts. In general, there are two approaches in modeling the deformation of flexible parts: Geometrical: In this approach, single or multiple control points or shape parameters are used for manual adjustment of a flexible part in order to apply changes in its shape and model its desired configuration. Generally, Geometrical methods rely on the designer’s skill rather than on physical principles, are computationally efficient, and aim to find a possible deformation in a reasonable time although they often produce results not matching with desired configurations. A common method for deformation of geometrical properties is the Functional method, in which a shape is expressed as a Bezier, B-spline, rational B-spline, or NURBS curve, and the designer reshapes the curve by adding, deleting, moving, or changing weights of the control points [3]. A much more powerful method is the Free-Form Deformation (FFD), which deforms the shape of a flexible object by changing the space that encompasses it. This is done by applying translation, rotation, bending, and other transformation matrices to change the model of the object. More complex deformations can be derived via a combination of these matrices [30]. Physical: In this approach that has been mostly used by animators and graphic designers, some sort of integration mechanisms and physical principles (e.g., material characteristics, environmental constraints, and external applied forces on a part) are used to compute the shapes or motions of flexible objects. A flexible part is shown as an entity with many degrees of freedom (DOFs), one for each deformable point on the surface of the object. The main disadvantage of Physical methods is the increase of DOFs to an uncontrollable level; thus, a balance must be made between the required time and the degree of realism. The most commonly used strategies for building deformable objects with physical properties are: (1) Mass-spring; (2) Finite Element Methods; and (3) PointBased systems [12]. Assembly sequence planning of flexible parts has been rarely addressed in the literature. Wolter and Kroll [40] discussed different types of operations on particularly string-like parts such as threads, ropes, wires, cables, and hoses. They described the state of an assembly by a set of relations among the features of its component parts together with the basic algorithms to determine the possibility of performing a particular operation on string-like assemblies. However, that work is not applicable to 2D and 3D flexible parts, and no method for ASP is proposed. A few other works deal with creating contact states of the flexible parts of an assembly. In that approach, information for the contact states is extracted from the geometrical models of the parts and the assembly and is used to construct a graph representation proper for assembly sequence and path planning. In the Contact State Graph, each node represents a valid contact state between two parts, and each arc between two nodes implies the adjacency of their corresponding contact states. Based on this graph, relative assembly paths can be planned as the sequences of contact state transitions connecting the starting state to the goal state. Remde et al. [28] identified different possible contact states between a linear deformable object and a rigid polyhedral body, and listed the feasible transitions between these states. A further elaboration on this formalism characterizes contact states by their stability and defining contact state transition classes [1]. Also in [21] the contact forces of an elastic tube during assembly operations are modeled based on different types

129

of contact states with the effects of deformation and friction taken into account. None of the above works, however, directly focus on ASP with flexible parts, and as a result this field of research has remained largely untouched. This paper aims to fill in this gap, and proposes a new formal framework for defining assemblability of flexible parts, in which concepts like Assembly Interference Matrix (AIM) and Assembly stress matrix (ASM) for flexible parts are mathematically defined from not only geometrical, but also physical and mechanical aspects, as presented and exemplified in Section 2. Another contribution of this work is solving the ASP problem for products with rigid and flexible parts by customizing the components of the Scatter Search (SS) metaheuristic algorithm, tuning its parameters to their best values through a systematic combination of TOPSIS and Taguchi methods as described in detail in Section 3. In Section 4 the presented method is extensively compared with GA, MA, SA, Breakout Local Search (BLS), Multistart Local Search (MLS), and Iterated Local Search (ILS) methods by solving three sample assemblies with different geometries, dimensions, and complexity via the detailed statistical analyses. In this section, a computational and statistical analysis of flexible vs. rigid modeling is also presented. Finally, concluding remarks are presented in Section 5.

2. Assembling flexible parts Usually, assembling a flexible part requires exerting forces to push it through a narrow passage or establish a latching contact with a mating part. Such a force, however, must strain the part within its elastic limit to avoid permanent deformation (and probably damage) of the part. In order to consider the deformability of flexible parts in assembly sequence planning, their elastic behavior under variable forces must be studied first. In this section, new formal definitions and mathematical symbols for showing interference relations between parts of an assembly and the amount of compressive stress needed for assembling flexible parts are presented, and then the process of creating the Assembly stress matrix (ASM) – which is a generalization of the Assembly Interference Matrix (AIM) – using the Abaqus software is described step-bystep to sufficient details. Lastly, the ASM is constructed for a sample assembly with two rigid and four flexible parts. In this paper, deformability of the flexible parts of assemblies during the assembly process are simulated by the AbaqusTM software, which is the reference software for Finite Element Method (FEM) analysis, modeling of deformation behavior of mechanical parts, and Computer Aided Engineering (CAE). For the FEM analysis, the user must input the part’s geometrical dimensions, material, density, Young’s modulus, coefficient of friction, types of elements used for modeling its deformation behavior, as well as the direction of exerted forces, into the Abaqus software. The software then simulates the part’s motion along the given force direction and calculates the exerted stresses and deformations for all finite elements of all contacting parts, which are used for constructing the Assembly stress matrix. It is noted that part variations which affect the assemblies are not considered in this paper.

2.1. Assembly stress matrix (ASM) In order to construct the Assembly stress matrix, we precalculate all possible collisions between any two parts at all directions and organize this information as an n × n × 2 m matrix called Assembly Interference Matrix (AIM). Each entry of the AIM k , which takes value 0 if part p is denoted by a binary variable Ii,j j does not block the movement of part pi along direction dk , and takes 1 otherwise. This variable is formally defined in (1) based on

130

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Table 1 Symbols used for denoting the geometric characteristics of parts during their assembly. Symbol

Description

BB(i )

The bounding box of the parts j (j = 1, . . ., i − 1) assembled prior to i which has axes parallel to the main Cartesian directions x, y, and z The Assembly Vector of part i along direction dk , starting from a point on the boundary of BB(i ) and ending at the final configuration of i in the assembly A The closure of the volume occupied by part i in its final configuration The swept volume of part i along AV(i , k) in the workspace

AV(i , dk )

cl(V(i )) SVAV (i ,dk ) i

p6 p3

p5

p1

k Ii,j =

0

AV (pi ,k)

if SVi

∩ cl(V (pj )) = ∅

,

Fig. 1. A puzzle-like assembly with four flexible (p1 , p2 , p5 , p6 ) and two rigid (p3 , p4 ) parts.

of the part. Accordingly, each element of the ASM (in megapascals, MPa) is calculated as:

(1)

1 otherwise k Si,j =

k we only need to check whether It is noted that for calculating Ii,j pj is on the way of pi when moving from its initial position toward its final position in the assembly along direction dk . Therefore, existence of other parts along this path is not considered and so there is only a unique AIM matrix for an assembly. For all pairs of parts pi and pj (i = 1, . . ., n; j = 1, . . ., n) and for all directions dk ∈{−x, +x, −y, +y, −z, +z} in 3D, the AIM is an n × (n × 2 m) matrix as defined in (2). Now we can define a solution s = (1 , d1 ), (2 , d2 ), . . ., (n , dn ) k = 0, meaning that there to an ASP problem is feasible if for all j < i, Ii,j is no already-assembled part j blocking i along its motion along direction di . The ASM is similar to the AIM both in dimension and k in the AIM are replaced with S k in the structure, except that all Ii,j i,j ASM. In constructing the ASM, the results of FEM analysis for the flexible parts performed in the Abaqus play an important role. p1

 ⎡

−x

+x

−y

0

0

0





+y

−z

+z

0

0

0



...

⎧ 0 ⎪ ⎪ ⎨

k + k ˆ i,j ˆ i,j

⎪ ⎪ ⎩ m

+x In, 1

−y

In, 1

+y

In, 1

−z In, 1

+z In, 1

...

if

k = 1) ∧ ( k ≤ Y ) ∧ ( k ≤Y) (Ii,j ˆ i,j ˆ i,j i j

if

k (Ii,j

=

k 1) ∧ ((ˆ i,j

>

k Yi ) ∨ (ˆ i,j

,

(3)

> Yj ))

In this subsection it will be shown how the ASM of a puzzle-like product with four flexible and two rigid parts (Fig. 1) is constructed. The materials and mechanical properties of the parts are shown in Table 3. First the AIM of the product is created by considering pure geometrical blockages, as if all parts were rigid. This is done by

pn

+y

+x

−y

−x I1, n

+x I1, n

I1, n

I1, n

−z I1, n

+z I1, n

0

0

0

0

0

0

+y

−z



−x

Suppose that for FEM analysis, part pi is subdivided into li polygonal or polyhedral small elements (for 2D and 3D parts, respectively), and ept i is the t-th element of part pi . Let us define some symbols with regard to the stresses exerted on parts during their assembly, as shown in Table 2. In order for a flexible part pi to remain within its elastic limit, k ≤ Y (∀j) must hold, in which Y is the yield stress the relation ˆ i,j i i

k = 0) (Ii,j

2.2. ASM construction



−y

if

in which  m is the maximum stress that can be exerted by the k is calculated for all pairs of parts p and p (i = 1, . . ., assembler. Si,j i j n; j = 1, . . ., n) and for all directions dk ∈ {−x, +x, −y, +y, −z, +z}, and thus the ASM is constructed.

+z

⎤

⎢ I −x I +x I −y I +y I −z I +z . . . I −x I +x I −y I +y I −z I +z ⎥ ⎢ 2, n 2, n 2, n 2, n ⎥ 2, n 2, n AIM = ⎢ 2,1 2, 1 2, 1 2, 1 2, 1 2, 1 ⎥ ⎣ ⎦ ... ... ... −x In, 1

+y

p2

+x

geometric characteristics of the parts during their assembly shown with symbols presented in Table 1.



p4

(2)

arranging  the coordinates of the vertices of assembly parts in an m × vi matrix. Then, each row of the AIM is constructed by translating its corresponding part (say pi ) from a point on the boundary of the (axis-aligned) bounding box of each of the other n − 1 parts to its final position in each of the 2 m directions. Translation along a certain direction (say dk ) is realized by step-wise movements with small increments, and each new position of the part is checked for

Table 2 Symbols used for denoting the exerted stresses on parts during their assembly. Symbol

Description

k i,j (t)

The stress applied on element ept during the assembly of part pi along direction dk after the part pj has been assembled

k i,j (t) k k ˆ i,j = max{i,j (t), k k ˆ i,j = max{i,j (t), k k ˆ j,i = max{j,i (t), k k ˆ j,i = max{j,i (t),

The stress applied on element ept during the assembly of part pi along direction dk after the part pj has been assembled

i

j

t = 1, . . ., li }

The maximum stress exerted on part pi during the process of its assembly along direction dk after part pj has been assembled

t = 1, . . ., lj }

The maximum stress exerted on part pj during the process of assembling part pi along direction dk after part pj has been assembled

t = 1, . . ., lj }

The maximum stress exerted on part pj during the process of its assembly along direction dk after part pi has been assembled

t = 1, . . ., li }

The maximum stress exerted on part pi during the process of assembling part pj along direction dk after part pi has been assembled

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

131

Table 3 Mechanical properties of the rigid and flexible parts of the assembly shown in Fig. 1. Type

ID

Material

Density (kg/m3 )

Yield stress, Y (MPa)

Young’s modulus, E (GPa)

Coefficient of friction

Finite element typea

Rigid

p3 p4

Aluminum Stainless Steel

2700 8190

240 520

68.9 180.0

0.10 0.10

C3D8R

Flexible

p1 , p2 p5 , p6

Nylon6 PVC (Polyvinyl Chloride)

1130 1300

45 52

3.0 3.3

0.05 0.05

C3D4

a

The finite element type of rigid parts is C3D8R (8-node linear brick element), and the flexible parts are of type C3D4 (4-node linear tetrahedron element).

collisions with other parts. In case of any collision with any other part (say pj ), the element in row pi and column dk of part pj is set to 1. This process is repeated for all parts of the product. For the puzzle-like assembly, construction of the AIM (shown below) took 37 s.



p1 p2 p AIM = 3 p4 p5 p6

p1



 

−x + x − y + y ⎡ 0 0 0 0 ⎢ 1 1 1 1 ⎢ 1 0 0 1 ⎢ ⎢ 0 0 0 1 ⎣ 0 0 0 1 0 0 0 1

p2



 

−x + x − y + y 1 1 1 1 0 0 0 0 1 0 1 1 1 0 0 1 1 0 1 1 1 0 0 1

p3



 

−x + x − y + y 0 1 1 0 0 1 1 1 0 0 0 0 0 1 1 1 0 1 0 1 0 1 0 1

support planes or fixtures are used. All forces and boundary conditions are created and exerted in Load module of the software, which has an option for setting their values not to exceed the parts’ elastic limits. p4



 

−x + x − y + y 0 0 1 0 0 1 1 0 1 0 1 1 0 0 0 0 0 1 1 0 0 1 0 1

p5



 

−x + x − y + y 0 0 1 0 0 1 1 1 1 0 1 0 1 1 0 1 0 0 0 0 0 1 1 1

p6





−x + x − y + y ⎤ 0 0 1 0 0 1 1 0 ⎥ 1 0 1 0 ⎥ ⎥ 1 0 1 0 ⎥ ⎦ 1 0 1 1 0 0 0 0

Now the ASM matrix is constructed as follows: for any combination of pi , pj , and dk in AIM for which Ii,k j = 0, the corresponding k in ASM is set to 0, but for those p , p , and d for which element Si,j i j k k = 1 and one of (or both) parts are flexible, their corresponding Ii,j assembly operation must be simulated in the Abaqus software to check if pj blocks pi along dk . Such a simulation is performed via different modules of the software, as follows:

1. The two parts are created in Part module of the software as deformable or rigid objects. 2. Density and mechanical properties of the parts (e.g., Young’s Modulus, Poisson’s Ratio, and Yield Stress) are assigned to the parts in Property module. 3. Part pj is placed at its fixed position in the final assembly, and part pi is placed at a point on the outer boundary of the bounding box of part pj which coincides with the line with direction dk along which part pi is supposed to translate toward its final position in the assembly. This is done in the Assembly module of the software. 4. In Step module, two types of steps must be defined: an ‘Initial step’ in which static and general conditions of the parts are set; and a ‘Dynamic step’ in which dynamic forces and motions are exerted and applied. 5. In Interaction module, all surface–surface contacts and interactions between parts are added to the Initial step. These contacts, however, may change in subsequent Dynamic steps. 6. Three forces are exerted on part pi in the Dynamic step: (1) the body force f1 exerted by the assembler along dk ; (2) the friction force f2 applied along contact surfaces with part pj ; and (3) the gravity force f3 applied to all parts along direction −y. The friction and gravity forces are applied to the part pj as well, and since it must be fixed (unlike part pi that translates under the body force), its displacements and rotations are constrained by the friction and contact forces exerted by all of its adjacent parts (that is, we presume that they are assembled before pj and thus prevent its free movements). If all previously-assembled parts do not suffice to constrain all motions of pj (like when it lies on the outer boundary of the final assembly), then auxiliary virtual

7. All parts are seeded and meshed in Mesh module. The shape of elements of flexible parts is hexahedron and the ‘sweep technique’ is used for meshing those parts. On the other hand, the shape of elements and meshing technique for rigid parts are tetrahedron and ‘free’, respectively. 8. In the Job module of the software, the type of the created job is set to ‘full analysis’. Other job types are ‘recover’ and ‘restart’, which we did not use. 9. The results of the FEM analysis is generated and reported by the software after a certain amount of time.

+x As an example, for calculating S1,2 for the flexible parts p1 and p2 in Fig. 1, the translation of p1 along the direction +x is simulated in Abaqus by exerting forces until it encounters the part p2 which is assumed to be fixed at its final assembly position. At the moment of encounter, the vector sum of the forces applied to the parts causes parts p1 and p2 to go under stress, which is calculated for each finite +x element of the two parts. The value of S1,2 is then set to the sum of the largest stresses in all elements of each part, which was 53 MPa for parts p1 and p2 . Fig. 2 shows simulation snapshots of assembling the part p1 along the direction +x after that the part p2 has been assembled. The magnitudes of applied stresses on different elements of the parts (the values of stresses of each contour color) are represented by the color scale (according to the contour plot legend) shown in Fig. 2(a), which ranges from dark blue, light blue, green, yellow, orange, to red for indicating the lowest to the highest stresses. The translational motion of part p1 continues until either it reaches its desired assembly configuration while the condition k ≤ Y ) ∧ ( k ≤ Y ) remains satisfied, or the parts undergo ˆ i,j (ˆ i,j i j

plastic deformation before the part p1 reaches its desired assembly configuration. The ASM of the studied assembly is constructed using the results of the FEM analysis as follows:

132

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Fig. 2. Simulation snapshots of assembling the flexible part p1 after part p2 along the direction +x: (a) the initial position of the part p1 (on the left); (b)–(d) intermediate positions of the part p1 . The colors of part elements indicate the accumulative exerted forces according to the color legend. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

p1

p1 p2 p ASM = 3 p4 p5 p6

  −x + x − y ⎡ 0

0

0

⎢ 49 m m ⎢ m 0 0 ⎢ ⎢0 0 0 ⎣ 0 0

0 0

0 0

 

p2



 

p3



+ y −x + x − y + y −x + x 0 m 53 m m 0 m m 0 0 0 0 0 m m m 0 m m 0 0 m m 0 0 m 0 m m m 0 m m 0 m m m 0 0 m 0 m

−y m m 0 m 0 0

 

+y 0 m 0 m m m

3. The assembly sequence planning algorithm The goal of Assembly sequence planning is finding the best sequence for assembling the rigid and flexible parts of an assembled product while minimizing one or more performance criterion(s) such as number of assembly direction changes, maximum required stress for accomplishing the assembly operation, and number of satisfied geometrical (precedence) constraints. As mentioned earlier, the ASP is an NP-hard optimization problem and thus exact mathematical programming approaches cannot solve large-scale instances of the problem within reasonable times. Therefore, metaheuristic algorithms have found wide applications in this realm. On the other hand, to the best of our knowledge, all previous works on ASP solution have addressed assemblies with rigid parts only, and so efficient solution of ASP with merely flexible or a combination of flexible and rigid parts is an open research theme. The present paper, however, addresses this new problem possibly for the first time by implementing a metaheuristic approach. It is noted that our method is designed for solving ASP of solid (2D and 3D) flexible parts, and not string-like parts such as wires and cables.

−x 0 0 m 0 0 0

p4



 

p5







p6





+ x − y + y −x + x − y + y −x + x − y + y ⎤ 0 m 0 0 0 m 0 0 0 m 0 m m 0 0 m m m 0 m m 0 ⎥ 0 m m m 0 m 0 m 0 m 0 ⎥ ⎥ 0 0 0 m m 0 m m 0 m 0 ⎥ ⎦ m m 0 0 0 0 0 m 0 m 69 m 0 m 0 m 62 m 0 0 0 0

Generally, in designing metaheuristic algorithms for solving optimization problems, two types of search strategies may be adopted: exploitation (intensification), and exploration (diversification) [33]. Algorithms emphasizing the first strategy are powerful in intensifying the search locally around local optima and are called Single-solution-based metaheuristics (or S-metaheuristics) since a single solution is iteratively improved until a satisfactory solution is obtained. The second class of algorithms introduce diversity in exploring the search space by simultaneously and iteratively improving multiple solutions and finally select the best solution among them, which are called Population-based metaheuristics (or P-metaheuristics). Local Search, SA, and Tabu Search (TS) are examples of S-metaheuristics, and GA, PSO, ACO, AIS, and HS are examples of P-metaheuristics. There are also hybrid algorithms, which usually nest an S-metaheuristic inside a P-metaheuristic, such as GA-SA and GA-TS. In this paper, Scatter Search (SS) algorithm is selected for solving the ASP problem. Scatter Search is an evolutionary and population-based metaheuristic first introduced in [13]. The algorithm starts with generating an initial population by considering

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Creating an initial population

133

Updating the Reference Set

Improving the initial population

Selecting the Reference Set

Generating Subsets

Improving the obtained solutions

Solution recombination

Fig. 3. Illustration of the search mechanism in the Scatter Search algorithm.

two key factors: quality and diversity of solutions. Several strategies can be applied to get a population with such properties. In general, greedy procedures are used to diversify the search, and a simple heuristic procedure is applied to the solutions to refine them. Next, some high quality and diverse solutions are selected to form a ‘Reference Set’, the elements of which are then grouped in (usually two-member) subsets and recombined to generate new solutions. Afterwards, these solutions undergo improvement through a selected neighborhood-searching procedure and integrate with the current Reference Set in order to update it toward better quality solutions. The above process iterates until a stopping criterion is satisfied. The final best solution is then selected from among the last population of improved solutions. Fig. 3 illustrates the main components of the Scatter Search algorithm. The reasons for selecting the SS algorithm for our ASP problem are: (1) it benefits from both the diversity inherent in generating populations of solutions, and the quality achieved by improving a single solution through local search, thus securing proper exploration and exploitation of the search space; (2) it implements solution improvement procedures on both the initial population of solutions, and the evolving populations throughout the search process. The following subsections describe details of our implementation of the Scatter Search algorithm for solving the ASP problem with flexible and rigid parts.

by a quality metric. In fact, these two properties guide us in selecting the ‘best’ neighbor(s) of a specific solution during the search process and hence converging to desired final solutions. Since all solutions at each iteration of the Scatter Search need to be checked for feasibility, we define the feasibility of a solution according to the definition of the ASM in (3). A solution s = (1 , d1 ), (2 , d2 ), . . ., (n , dn ) is feasible if either there is no already-assembled part j blocking i along di , or the stress required for performing this assembly is less than the maximum stress exerted by the assembler, expressed as:

3.1. Solution representation

MAS(i ) =

∀i = 2, . . ., n,

d i S < m , i ,j

j = 1, . . ., i − 1.

(4)

In this paper, the objective function is to simultaneously minimize the total number of Assembly Direction Changes (ADC) and the amount of Maximum Applied Stress (MAS) for performing the assembly operation. ADC(i ) is a binary variable indicating the change of assembly direction between parts i−1 and i in the solution string:



ADC(i ) =

0

if di = di−1

(5)

1 otherwise

Also, MAS(i ) is computable using the solution vector and the ASM as explained in Section 2.2: i−1 

di S i ,j

(6)

j=1

The first step in designing a metaheuristic is proper representation (encoding) of solutions. A solution s to the ASP problem can be represented by the assembling sequence and directions of the n parts of the assembly A in the form of a row vector of ordered pairs (i , di ), in which i ∈{p1 , p2 , . . ., pn } is the i-th element of the sequence and di ∈{−x, +x, −y, +y, −z, +z} represents the assembly direction of part i (for 3D assemblies). Fig. 4 shows a typical encoding of a solution as (p5 , −x), (p3 , +y), . . ., (p6 , +x), which indicates that part 1 (=p5 ) must move at direction −x toward the assembly, then part 2 (=p3 ) must move at direction +y toward the assembly, and so on. Since the total number of possible sequences of n parts is n! and the number of possible directions to assemble a particular part is 2m, then the size of the whole search space (including feasible and infeasible solutions) equals to n!(2m)n . Also, since the solution encoding has 2n cells, the maximum distance between two solutions in the search space (i.e., the diameter of the search space) is 2n, as all cells in two solutions may be different. 3.2. Fitness of solutions Solving an ASP problem requires finding assembly sequences and directions that are both feasible and high quality as measured

5

–x

3

+y

2

+y

1

–y

11

–y

7

Now the total fitness function can be defined as the weighted sum of the above criteria expressed as: Minimize

f (s) = w1

n 

ADC(i ) + w2

i=2

n 

MAS(i ).

(7)

i=1

Minimizing the first term in (7) is equivalent to minimizing the total number of assembly direction changes, and minimizing the second term is equivalent to minimizing the maximum applied stress to the assembly parts to perform the assembly operation. In order to bias the search toward feasible solutions, we set the values of the weighting factors as w2 = k × w1 . The value of k is adjusted later in Section 3.8. 3.3. Generation, diversification and improvement of initial solutions The Scatter search starts with a set of initial solutions, which are required to be as scattered as possible within the solution space. Thus, the use of a systematic procedure to generate such solutions is essential. Regarding the ASP as a permutation problem, the described method in [14] for generating diverse solutions is considered and modified, which operates as follows:

–y

9

+x

10

–x

8

Fig. 4. A typical encoding for assembly sequence of 11 parts in 2D.

+x

4

–x

6

+x

134

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

+y

1

2

–y

3

–y

4

+x

5

–x

6

+x

f = 4w1 + (7σm)w2

5

–x

2

–y

3

–y

4

+x

6

+x

f = 3w1 + (8σm)w2

Insertion position

1

+y

Fig. 5. Applying the Random Insertion operator to the pair (5 , d5 ).

1

+y

2

–y

3

–y

4

+x

5

–x

6

+x

f = 4w1 + (7σm)w2

2

–y

6

+x

3

–y

4

+x

5

−x

f = 5w1 + (6σm)w2

Insertion position

1

+y

Fig. 6. Applying the Effective Insertion operator to the pair (6 , d6 ).

A solution s = (1 , d1 ), (2 , d2 ), . . ., (n , dn ) is used as a starting point (seed) to generate n − 1 subsequent permutations. Let us define the sub-sequence s(h: q) where q is a positive integer between 1 and h, such that: s(h : q) = (q , dq+rh ), (q+h , dq+(r−1)h ), (q+2h , dq+(r−2)h ), . . ., (q+rh , dq )

(8)

where r is the largest nonnegative integer such that q + rh ≤ n. Finally, the solution s(h) for 2 ≤ h ≤ n is defined as: s(h) = s(h : h), s(h : h − 1), . . ., s(h : 1)

(9)

For example, assume h = 3 and a seed permutation in the ASP problem is s = (1, +x), (2, −y), (3, −y), (4, −x), (5, −y), (6, +y), (7, +y), (8, −x). The recursive application of s(3: q) for q = 3, 2, 1 results in subsequences s(3:3) = (3, +y), (6, −y), s(3:2) = (2, −x), (5, −y), (8, −y), and s(3:1) = (1, +y), (4, −x), (7, +x). Hence, s(3) = (3, +y), (6, −y), (2, −x), (5, −y), (8, −y), (1, +y), (4, −x), (7, +x). About 20% of the initial Population Size (PS) is generated by the above method and the rest by random permutation. After generating diverse initial solutions, their quality is improved by applying a new proposed insertion operator called Effective Insertion, which is a variation of the Random Insertion operator and aims to produce solutions with smaller total Maximum Applied Stress. The Random Insertion operator randomly selects a pair in the assembly sequence and relocates to another random position in the string, such that all pairs after that position shift one place to the right, as shown in Fig. 5. The neighborhood size for any solution through applying this operator is n × n = n2 , since n pairs can be selected as the start, and n positions (including the position after the last pair) can be selected as the target of the relocation. The Random Insertion may or may not decrease the total Maximum Applied Stress. So, to make it more goal-directed, the Effective Insertion operator is proposed which acts more ‘intelligently’ since it selects the insertion position by also considering the Maximum Applied Stress caused by the inserted part rather than by just choosing it randomly. Each application of the Effective Insertion operator is as follows: for i = 2, . . ., n calculate MAS(i ), find the part with the maximum MAS value, and insert it before the first part in the sequence that blocks its assembly operation. The Effective Insertion operator is applied only to infeasible solutions, and considering its caused relatively small perturbation in the solution, we apply this operator a few times until either the (infeasible) solution becomes feasible or an upper limit is reached. For example, for the assembly in Fig. 1 and the solution (1, +y), (2, −y), (3, −y), (4, +x), (5, −x), (6, +x), the Maximum Applied Stresses for all six parts are MAS(i ) = [0, m ,

 m ,  m ,  m , 3 m ] and so the part 6 with the maximum MAS value is selected and inserted before the first part in the sequence that blocks its assembly operation, which is part 3 , as shown in Fig. 6. Accordingly, for the obtained solution, MAS(i ) = [0, m , 0, 2 m ,  m , 2 m ] with a total Maximum Applied Stress of 6 m which is smaller than the previous 7 m . It is worthy to note that performing an Effective Insertion does not guarantee an improvement in the fitness function; however, our extensive experiments showed that it reduces the total Maximum Applied Stress far better than the Random Insertion operator. For feasible solutions or those that remain infeasible after applying multiple consecutive Effective Insertions, a number of Random Flip operators are applied by arbitrarily selecting a pair in the sequence and randomly changing its assembly direction. This operator is repeated until either the quality of the solution improves or an upper limit is reached. As a result, a neighboring solution is generated which is checked for improvement in the fitness function. If so, then the initial solution is replaced by this new solution. The above Effective Insertion and Random Flip operators are applied on all solutions of the initial population of size PS in order to improve their quality. 3.4. Reference set selecting and updating The Reference Set is an important component in the Scatter Search method, and contains some elite solutions that have sufficient diverse and high-quality solutions. It is suggested in the literature (e.g. in [20]) to collect the elements of the Reference Set by selecting the best and diverse solutions to avoid similarities in solutions. In our proposed method, the set of improved solutions obtained from the improvement method (Section 3.3) is sorted based on their objective function values, and the Reference Set is constructed from two types of solutions: 80% of solutions with better objective function values from the top of the list (called RefSet1), and 20% of the other unselected solutions with the highest Hamming distance to the solutions in RefSet1 (called RefSet2). The size of the Reference Set is equal to B, as calculated below, in which 0 < ˇ < 1 and PS indicates the initial Population Size: B = [ˇ × PS].

(10)

3.5. Subset generation In this step, a fraction (20%) of all possible two-element subsets of the Reference Set is randomly generated in the following manner: 2% with both elements taken from solutions in RefSet2 (which are characterized by good diversity), 4% with an element

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Parent A:

1

+y

5

–x

2

–y

3

–y

4

+x

6

+x

Parent B:

4

–y

2

–y

6

+y

5

+x

1

+y

3

–x

Child a:

2

–y

3

–y

6

+y

5

+x

1

+y

4

+x

Child b:

6

+y

5

+x

2

–y

3

–y

4

+x

1

+y

135

Fig. 9. Illustration of the Two-point Crossover operator. Fig. 7. An example of parents and children in the Partially Matched Crossover (PMX).

3.7. Improvement method

Parent A:

1

+y

5 –x

2

–y 3 –y 4

+x 6

+x

Parent B:

4

–y

2 –y

6

+y 5 + x 1

+y 3

–x

Child 1:

4

–y

5 –x

2

–y 3

+y 1 +y

–y 6

Fig. 8. Illustration of the Order Crossover (OX) operator.

taken from solutions in RefSet1 (which are characterized by good quality) and the other from solutions in RefSet2, and 14% with both elements taken from solutions in RefSet1. In this way, a mixture of high quality, high-diversity, and good quality-diversity pairs is created aiming to attain a proper intensification and diversification of the problem space. These pairs will be further manipulated in the next step for creating new and hopefully better solutions.

3.6. Solutions recombination After producing a subset of solutions by the subset generation method, they are recombined to create new solutions using different combination methods, as discussed in [23]. In this paper, the following three procedures are considered as recombination operators: Partially Matched Crossover (PMX): In this operator, in general, two genotypes (solution encodings in the form of strings of n (part, direction) pairs) are selected as parents A and B, and two crossover positions are picked randomly along the solutions. Then, all the chromosomes (elements) of Parent A lying between these two points are exchanged with the chromosomes of Parent B at the same absolute positions, and vice versa. The crossover is finished by placing the remaining elements into the vacant positions of the offspring from left to right, in order of their appearance in the corresponding Parent. For example, for the ASP strings in Fig. 7, taking the Parents A and B, the two crossover limits are fixed at 2nd and 5th positions, and the pairs between these positions must undergo exchange. Ordered Crossover (OX): In this method only one offspring is generated from two parents. First a substring from the Parent A is selected randomly and an offspring is produced by copying the substring into its corresponding position. Then these selected elements are deleted from the Parent B. The resulting sequence contains the elements that the offspring needs. The crossover is finished by placing the remaining elements into the vacant positions of the offspring from left to right, in order of their appearance in the Parent B. The procedure is demonstrated in Fig. 8. Two-point Crossover: In this operator, two points are randomly selected, and the elements outside the selected two points are inherited from parent A of the child, and the other elements are placed in the order of their appearance in the parent B, as shown in Fig. 9. In this method (similar to the OX crossover) only one offspring is generated from two parents.

In order to improve each candidate solution in the set of new solutions (children), the Effective Insertion and Flip operator are applied to it, as described in Section 3.3. Afterwards, a pool of N solutions is obtained, that is the collection of parents (those in the Reference Set) and improved children. The next Reference Set will be updated by selecting high quality and diverse solutions from the above pool with the same methodology explained in Section 3.4, while its size B is maintained. The Updating-Recombination-Improvement procedure (the cycle in Fig. 3) is repeated until the stopping criterion is satisfied, which is when the number of iterations reaches an upper bound like Itermax , which we set to 30. 3.8. Parameter tuning As stated before, parameters of metaheuristic algorithms have significant effects on their efficiency and effectiveness in solving a particular problem. As there may be many possible values for parameters of a given algorithm, using an appropriate parameter tuning approach to find the best choice from many alternatives can produce better solutions for the problem. In off-line parameter tuning (that is, before the execution of the algorithm), the designer tunes one parameter at a time and determine its optimal value empirically. As a result, interactions between parameters are not considered and so finding the optimal setting (even through an exact optimization method) is not guaranteed. Therefore, the Experimental Design (or Design of Experiments (DOE)) approach in general, and the Taguchi method in particular, is used in this paper to investigate different scenarios of parameter settings and to select the best values of parameters concurrently. Here we seek to tune the parameters of the SS algorithm to optimize a combination of the following criteria, in which N is the total number of solutions in the last pool of the run experiment, and ik is the i-th part in the sequence of the k-th solution of the last pool: (1) minimizing the total number of Fitness Function Evaluations (FFE) during the whole search, with a weight of 0.15, ¯ = (2) minimizing the average total Maximum Applied Stress, MAS

N n i−1 k=1

i=2

S

d k  i

j=1 k ,k i

/N, with a weight of 0.30,

j

¯ = (3) minimizing the average Assembly Direction Changes, ADC N n k ADC(i )/N, with a weight of 0.10, k=1 i=2 (4) minimizing the objective function f(s* ), with a weight of 0.45.

In order to tackle these four objectives (which have different dimensions) simultaneously, the Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) method is used, which is a multi-criteria decision analysis method originally developed by [8] to assess alternatives prior to multiple-attribute decision making. In the TOPSIS, all the criteria are normalized to equalize their dimensions and the prespecified weights of the criteria are applied. Then the ideal solution V+ (the one which has the best level for all attributes (parameters)) and negative-ideal solution V− (the one

136

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

by considering number of levels − 1 = 3 degrees of freedom for the first three parameters, 2 degrees of freedom for the fourth parameter, and 1 degree of freedom for the total mean, summing to at least (3 × 3) + (2 × 1) + 1 = 12 run experiments. Therefore, L 16 (43 × 3) is chosen as the orthogonal array [26]. For choosing the best values of parameters, the average TOPSIS index (C¯ i ) for each parameter must be calculated and the level with the highest C¯ i value (i.e., closest to 1) among them is selected as the best level for that parameter. The result of parameter tuning is shown in Table 5, which indicates that the best levels for the PS, k, ˇ, and recombination type are 3, 3, 3 and 3, respectively (shown in bold), meaning that PS = 150, k = 9, ˇ = 0.15, and Recombination type = Two-point.

Table 4 Parameter levels in the Scatter Search algorithm. Parameter

Index of levels

Levels

PS (the initial Population Size)

1 2 3 4

50 100 150 200

k (=w2 /w1 in the objective function)

1 2 3 4

3 6 9 12

ˇ (ratio of the RefSet size to PS)

1 2 3 4

0.05 0.10 0.15 0.20

Recombination type for generating new solutions in the RefSet

1 2 3

PMX OX Two-point

4. Experimental results In this section, the results of implementing the Scatter Search algorithm for ASP of rigid and flexible parts are presented. The program was coded in MATLAB and run on a PC with an Intel Core i7® 3.2 GHz CPU and 16 GB of RAM. Since this is the first time flexible parts are considered in the ASP problem, there were no benchmark problems or methods in the literature that could be used to assess and evaluate the effectiveness and efficiency of our developed method. As a result, we selected three typical assemblies with various complexities and rigid and flexible parts, and solved them by the Scatter Search and numerous other metaheuristic algorithms. The selected assemblies and their specifications are as follows:

which has the worst attribute values) according to each alternative are computed. The ‘separation’ of each alternative from the ideal and negative-ideal solutions (S+ and S− , respectively) are calculated, based on which the TOPSIS index Ci = Si − /(Si + + Si − ) is calculated. The closer Ci is to 1, the higher is the priority of the i-th run experiment. Because of the multiple levels of each parameter, C¯ i is calculated as the mean of relative closeness to the ideal solution for each parameter per level. Finally, the best parameter setting is selected as the one with largest C¯ i . Details of the TOPSIS methodology are explained in [35] and [24]. Influential parameters in the SS algorithm for the ASP problem include (i) the initial population size PS, (ii) the constant coefficient k that is the ratio of the weights in the objective function (7), (iii) the constant factor ˇ in (10), and (iv) the recombination type in Section 3.6. Possible levels of the above parameters are set as indicated in Table 4. Table 4 suggests that 4 × 4 × 4 × 3 = 192 experiments are required for the full factorial design (i.e., all possible experiments). However, considering the computational cost and time, and based on statistical theories, there is no need to test all combinations of the factors. Instead, the Taguchi method is used to design only meaningful and effective experiments. But first the number of degrees of freedom should be calculated in order to select an appropriate Taguchi orthogonal array. The number of degrees of freedom for selecting the appropriate Taguchi orthogonal array is obtained

1. A mockup of the Soviet medium tank T-34/76, with 54 rigid and 9 flexible parts [45]. The assembling process of parts is depicted in Sections 1–9 of Fig. 10, and the parts and their ID’s are described in the table below the figure, in which flexible parts are marked by an asterisk (*), 2. A snap fit electrical connector assembly (Fig. 11(a)) with 4 rigid and 4 flexible parts (p3 , p4 , p5 and p6 ) [2], 3. A compact box with snap fit connectors (Fig. 11(b)) with 10 rigid and 10 flexible parts (p1 to p10 ) (designed by us). In order to evaluate the efficiency and effectiveness of the developed Scatter Search algorithm, we solved the above assembly sequence planning problems by it and several other metaheuristic methods, including (1) Genetic Algorithm (GA) [11], (2) Memetic Algorithm (MA) [11], (3) Simulated Annealing (SA) [17], as well

Table 5 Taguchi Orthogonal arrays of L 16 (43 × 31 ). No.

Population Size

k

ˇ

Recombination type

Ci (TOPSIS index)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4

1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

1 2 3 4 2 1 4 3 3 4 2 1 4 3 1 2

1 2 3 2 3 1 1 2 3 1 2 3 2 3 3 1

0.39669 0.54437 0.62958 0.57219 0.76176 0.49198 0.60815 0.63892 0.80025 0.51695 0.60546 0.65653 0.55733 0.74912 0.69552 0.55622

C¯ 1 = 0.53571 C¯ 2 = 0.62520 C¯ 3 = 0.64480 C¯ 4 = 0.63955

C¯ 1 = 0.62901 C¯ 2 = 0.57560 C¯ 3 = 0.63468 C¯ 4 = 0.60596

C¯ 1 = 0.56018 C¯ 2 = 0.61695 C¯ 3 = 0.70447 C¯ 4 = 0.56365

C¯ 1 = 0.51400 C¯ 2 = 0.58365 C¯ 3 = 0.71546

C¯

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Fig. 10. Mockup of the Soviet medium tank T-34/76 [45].

137

138

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Fig. 11. (a) Snap fit electrical connector assembly [2]; (b) compact box with snap fit connectors.

Table 6 Parameters of the implemented GA, MA, and SA algorithms. Algorithm

Parameter Description

Value

GA, MA [11]

PS: Population Size wc : Crossover probability Crossover type

100 0.8

wm : Mutation probability Mutation type SA [17]

 0 : Initial Temperature Neighborhood generation approach r(k): The number of non-improving solutions at a particular iteration k The cooling schedule

Partially Matched Crossover (PMX) 0.1 Exchange operator 2 Swap + Exchange [2k/n(n − 1)]

 k =  0 /1 + ln(r(k))

as other versions of those algorithms with parameters tuned by ourselves, including (4) Memetic Algorithm (MA1), (5) Genetic Algorithm (GA1), (6) Simulated Annealing (SA1), and also the following metaheuristics implemented by ourselves: (7) Breakout Local Search (BLS), (8) Iterated Local Search (ILS), and (9) Multistart Local Search (MLS). The parameters of the GA, MA, and SA algorithms are listed in Table 6 as originally proposed by their developers. Also, the GA1 and MA1 algorithms used the same objective function (7)

and initial solution and neighborhood generation methods as in the SS method, but employed the mutation operator in each iteration for diversification. The MA used the guided local search which is based on reducing the number of assembly direction changes by arranging the components whose assembly directions are the same. Also, the MA1 used the steepest descent approach as its local search method. In the BLS algorithm, the jump magnitude varied between the (Mmin , Mmax ) interval. Both the MLS and ILS algorithms used the hill climbing local search and randomly generated initial solutions. However, the MLS simply starts from a random solution in each iteration after 50 non-improving solutions, while the ILS iterates by conducting a new local search from a solution that is a perturbation of the last local optimum via the Inversion operator. In order to make a fair comparison, the parameters of the GA1, MA1, SA1, BLS, ILS, and MLS algorithms were also tuned by the same TOPSIS-Taguchi method as in SS and the best values were selected for each algorithm (indicated in Table 7 with boldface). The following performance measures were considered for comparing the abovementioned algorithms:

i. Minimizing the fitness value of the overall best solution in 30 runs of each algorithm, i.e., fmin = min{f (si∗ )},, in which f(si * ) is the best-found fitness value in the i-th run of the algorithm, ii. Minimizing the fitness value of the poorest solution among the best-found solutions of 30 runs of each algorithm, i.e., fmax = max{f (si∗ )},

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146 Table 7 Results of the TOPSIS-Taguchi-based parameter tuning method for all implemented algorithms. Method Parameter

Levels

SS

PS (Population Size) k (constant factor) ˇ (constant coefficient) Recombination type

50, 100, 150, 200 3, 6, 9, 12 0.05, 0.10, 0.15, 0.20 PMX, OX, Two-point

GA1

PS (Population Size) Mutation probability Mutation type Crossover probability Crossover type

50, 100, 150, 200 0.001, 0.004, 0.007, 0.01 Swap, Exchange, Insertion, Inversion 0.3, 0.5, 0.7, 0.9 PMX, OX, Two-point

PS (Population Size) Mutation probability Mutation type Crossover probability Crossover type

50, 100, 150, 200 0.001, 0.004, 0.007, 0.01 Swap, Exchange, Insertion, Inversion 0.3, 0.5, 0.7, 0.9 PMX, OX, Two-point

Initial temperature Neighborhood generation approach The number of non-improving iterations at a particular temperature (Equilibrium State) ˛ (the constant factor of the geometric cooling schedule)

500, 1000, 1500, 2000 Swap + Exchange, Swap + Insertion, Swap + Inversion 25, 50, 75, 100

Minimum number of jump magnitude (Mmin ) Maximum number of jump magnitude (Mmax ) Neighborhood generation approach Maximum number of iterations with non-improving local optima

0.05n, 0.1n, 0.15n, 0.2n

MA1

SA1

BLS

ILS

MLS

0.6, 0.7, 0.8, 0.9

0.75n, 0.8n, 0.85n, 0.9n Swap + Exchange, Swap + Insertion, Swap + Inversion 0.45n, 0.5n, 0.55n, 0.6n

Neighborhood generation approach Perturbation approach Acceptance Criteria

Swap + Exchange, Swap + Insertion, Swap + Inversion Swap, Exchange, Insertion, Inversion Strong selection, Weak selection, Probabilistic acceptance criteria

Neighborhood generation approach Acceptance Criteria

Swap + Exchange, Swap + Insertion, Swap + Inversion Strong selection, Weak selection, Probabilistic acceptance criteria

iii. Minimizing the average fitness of the best-found solutions in 30 ∗ 30 runs of each algorithm, i.e., favg = f (si )/30, i=1 iv. Minimizing the standard deviation (STD) of the best-found solutions in 30 runs of each algorithm, v. Minimizing the average gap between the fitness of the bestfound solution in each run of the algorithm and the best-known 30 solution fbest , i.e., AvgGap = (f (si∗ ) − fbest )/(30fbest ), i=1 vi. Maximizing the percent of feasible solutions (or success rate, SR) in 30 runs of each algorithm. This measure is important since usually more than half of the solutions in the search space of an ASP problem are infeasible, as they violate precedence relationships between parts of the assembly, vii. Minimizing the average runtime in 30 runs of each algorithm. For comparisons, we solved all the three benchmark problems introduced in Figs. 10 and 11 by the SS method and compared its performance measures to those obtained by the other nine algorithms, as reported in Tables 8–10, where best values in each column are indicated in bold. The stopping criterion of all algorithms was set to 70,000, 3000, and 30,000 total Fitness Function Evaluations (FFE) for the problems 1, 2, and 3, respectively. The best-known assembly sequence for the tank assembly (which is a feasible solution with 13 reorientations) was generated after 20

139

iterations of the SS method, with a better fitness value (=2596) than other algorithms. Also in the other two problems the SS outperformed other algorithms in all performance measures (except for runtime), thanks to the effective methods used for generating the initial solution, neighboring solutions, reference set selection and improved solutions. For example, Figs. 12 and 13 plot the average gap AvgGap and success rate SR of five out of ten algorithms that performed relatively well, which show the better performance of the SS method over the others.

4.1. Statistical analysis In order to statistically analyze and verify the significance of the outperformance of the SS algorithm, two nonparametric statistical tests have been conducted: (1) the “Friedman two-way ANOVA by ranks” test [10] for detecting any significant difference between behaviors of all algorithms, and (2) the “Kruskal–Wallis one-way ANOVA by ranks” test [19] for detecting a significant difference between the behaviors of each pair of algorithms. Note that parametric tests could not be used here because of their assumptions on independency, normality, and homoscedasticity of the variables [9]. The Friedman two-way ANOVA is a multiple comparisons nonparametric test and a counterpart to the parametric two-way ANOVA used to detect significant differences between the means of data from several groups. In our context, the Null Hypothesis (H0 ) of the Friedman’s test will be: “The ten algorithms SS, GA, GA1, MA, MA1, SA, SA1, BLS, ILS, and MLS have equivalent performances in solving the three benchmark ASP problems”, and the Alternative Hypothesis (H1 ) claims the opposite of H0 . An important result of the Friedman test is the relative performance ranks of the algorithms, the most superior one being ranked first (resp. last) for maximization (resp. minimization) criterion, and so on. The average relative rank j ri of each algorithm j in solving each benchmark problem i for 30 times is determined according to its performance in some selected criteria, and the global rank Rj of that algorithm is calculated by computing the average of its relative ranks. Another output is the p-values related to the similarity between algorithms, which indicate the maximum errors incurred when rejecting the H0 and accepting the H1 . For example, a p-value of zero indicates that the algorithms are definitely (with 100% of confidence) different. A commonly used confidence level for the p-value is ˛ = 0.05; thus, if the calculated p-value of a statistical test is larger than 0.05, we can reject the H1 regarding the difference of algorithms. Table 11 summarizes the numerical results of performing the Friedman test for the ten algorithms and the three benchmark problems. The p-values in the table show that performances of the algorithms in all criteria are significantly different at a confidence level of 100.0%. Also, values of the four important criteria vary from one problem to another. Now that there are significant differences between the performances of all algorithms, we conduct a Kruskal–Wallis test to characterize these differences for each pair of the algorithms. The Kruskal–Wallis one-way analysis of variance by ranks test is a pairwise comparisons nonparametric test and an analog of the parametric one-way ANOVA [19]. Here we conducted 9 × 4 × 3 = 108 tests using the best solutions of 30 runs of each algorithm as the test population of the respective method, and with a Null Hypothesis stating that “The performance of the SS algorithm is equivalent to the performance of algorithm i in problem j and criterion k”. If the Null hypothesis for a specific problem and a specific criterion is rejected (i.e., when p-value < ˛), then the SS algorithm stochastically dominates the compared algorithm. For each pair of SS and another algorithm, the Kruskal–Wallis test first sorts all the performance values regardless of their

140

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Table 8 The results of sequence planning of the Soviet medium tank T-34/76 mockup by SS, GA, GA1, MA, MA1, SA, SA1, BLS, ILS, and MLS methods. Algorithm

SS

GA

GA1

MA

MA1

SA

SA1

Best-found sequence f(s*)

(1, −z), (28, −z), (29, −z), (41, −z), (40, −z), (17, +x), (52, +x), (16, +x), (20, +x), (54, +x), (60, +x), (53, +x), (39, +x), (58, +x), (21, +x), (59, +x), (47, +x), (18, −x), (56, −x), (50, −x), (19, −x), (55, −x), (61, −x), (48, −x), (62, −x), (23, −x), (57, −x), (22, −x), (63, −x), (9, −x), (7, −x), (8, +x), (6, +x), (3, −y), (2, −y), (5, +y), (4, +y), (38, +x), (46, +x), (51, −x), (49, −x), (10, −z), (31, −z), (30, −z), (43, −y), (42, −y), (35, +y), (13, −z), (14, −z), (11, −z), (24, −z), (32, −z), (25, −z), (12, +y), (15, +y), (26, −z), (36, −z), (37, −z), (27, −z), (44, −z), (34, −z), (45, −z), (33, −z) (50, −z), (41,−z), (63, −z), (29, +y), (40, +y), (28, +x), (59, +z), (1, +z), (52, +z), (53, −x), (10, −y), (58, −z), (23, −z), (20, +z), (19, −x), (62 +z), (54, −z), (56, +y), (22, +x), (16, +x), (48, −x), (9, +y), (18, −z), (47, −z), (60, −z), (17, +x), (21, +x), (55, +z), (57, −x), (61, −x), (39, −y), (6, +x), (2, +x), (7, −x), (8, −x), (4, +y), (3, −y), (5, +y), (38, +x), (51, −x), (31, −z), (49, −x), (42, −y), (46, +x), (35, +y), (30, −z), (43, −y), (13, −z), (11, −z), (24, −z), (32, −z), (14, −z), (15, +y), (25, −z), (36, −z), (27, −z), (37, −z), (34, −z), (12, +y), (26, −z), (44, −z), (45, −z), (33, −z) (60, −y), (40, −y), (41, −y), (28, +z), (29, +z), (54, +x), (19, −x), (53, +x), (16, +x), (17, +x), (48, −x), (57, −x), (2, −y), (1, +x), (52, +x), (21 +x), (47, +x), (18, −x), (13, −z), (20, +x), (55, −x), (56, −x), (62, −x), (63, −x), (50, −x), (61, −x), (23, −x), (39, +x), (58, +x), (22, −x), (5, +y), (59, +x), (8, +x), (7, −x), (9, −x), (6, +x), (38, +x), (3, −y), (42, −y), (30, −z), (32, −z), (4, +y), (46, +x), (51, −x), (10, −z), (31, −z), (49, −x), (36, −z), (27, −z), (43, −y), (35, +y), (25, −z), (14, −z), (11, −z), (12, +y), (15, +y), (33, −z), (37, −z), (24, −z), (26, −z), (45, −z), (44, −z), (34, −z) (41, +x), (29, +x), (59, −x), (18, +x), (1, −x), (60, −z), (40, −z), (39, +x), (22, −x), (28, −z), (16, +x), (20, +y), (62, −y), (21, −x), (55, −z), (61, −x), (58, +x), (17, −x), (52, +y), (19, −y), (53, +x), (56, −x), (50, −z), (23, −x), (8, −x), (47, +x), (54, −x), (63, −x), (57, −x), (48, −x), (9, −x), (7, −x), (6, +x), (49, −x), (2, −y), (3, −y), (51, −x), (5, +y), (4, +y, (38, +x), (13, −z), (46, +x), (10, −z), (15, +y), (42, −y), (31, −z), (37, −z), (43, −y), (35, +y), (12, +y), (26, −z), (30, −z), (24, −z), (14, −z), (32, −z), (27, −z), (36, −z), (25, −z), (11, −z), (44, −z), (45, −z), (34, −z), (33, −z) (1, −z), (29, −z), (28, −z), (40, −z), (41, −z), (17, +x), (39, +x), (52, +x), (18, −x), (55, −x), (19, −x), (48, −x), (56, −x), (16, +x), (47, +x), (53, +x), (54, +x), (60, +x), (21, +x), (58, +x), (59, +x), (20, +x), (57, −x), (63, −x), (61, −x), (50, −x), (22, −x), (62, −x), (23, −x), (9, −x), (7, −x), (8, +x), (6, +x), (2, −y), (3, −y), (5, +y), (4, +y), (46, +x), (38, +x), (49, −x), (51, −x), (10, −z), (31, −z), (30, −z), (35, +y), (43, −y), (42, −y), (13, −z), (14, −z), (11, −z), (12, +y), (15, +y), (32, −z), (24, −z), (25, −z), (27, −z), (26, −z), (36, −z), (37, −z), (45, −z), (33, −z), (44, −z), (34, −z) (1, −z), (28, −z), (29, −z), (41, −z), (40, −z), (18, −x), (53, +x), (54, +x), (60, +x), (39, +x), (16, +x), (17, +x), (52, +x), (57, −x), (19, −x), (63, −x), (55, −x), (50, −x), (56, −x), (61, −x), (23, −x), (48, −x), (62, −x), (22, −x), (59, +x), (58, +x), (21, +x), (47, +x), (20, +x), (8, +x), (6, +x), (9, −x), (7, −x), (2, −y), (3, −y), (5, +y), (4, +y), (51, −x), (49, −x), (38, +x), (46, +x), (10, −z), (30, −z), (31, −z), (43, −y), (35, +y), (42, −y), (13, −z), (14, −z), (11, −z), (32, −z), (12, +y), (15, +y), (25, −z), (24, −z), (36, −z), (37, −z), (26, −z), (27, −z), (45, −z), (33, −z), (44, −z), (34, −z) (1, −z), (29, −z), (28, −z), (41, −z), (40, −z), (54, +x), (39, +x), (47, +x), (16, +x), (60, +x), (53, +x), (52, +x), (17, +x), (59, +x), (58, +x), (21, +x), (20, +x), (56, −x), (19, −x), (23, −x), (18, −x), (55, −x), (50, −x), (48, −x), (57, −x), (61, −x), (22, −x), (63, −x), (62, −x), (8, +x), (6, +x), (9, −x), (7, −x), (4, +y), (5, +y), (2, −y), (3, −y), (51, −x), (49, −x), (46, +x), (38, +x), (10, −z), (30, −z), (31, −z), (42, −y), (43, −y), (35, +y), (13, −z), (14, −z), (11, −z), (32, −z), (24, −z), (25, −z), (12, +y), (15, +y), (26, −z), (36, −z), (37, −z), (27, −z), (44, −z), (45, −z), (34, −z), (33, −z)

Fitness values

Efficiency metrics

fmin

fmax

favg

STD

AvgGap (%)

Success rate (%)

Runtime (s)

2596

2602

2600.60

1.84

0.18

100.0

139.79

224,957

431,623

303,227.10 69,066.55 11,580.55

0.0

507.77

109,488

228,360

204,088.40 42,964.46 7761.65

0.0

397.15

140,168

294,168

189,493.00 48,524.42 7199.42

0.0

309.26

2598

18,321

5324.13

5078.78

97.36

70.0

219.57

2599

2608

2603.20

3.05

0.28

100.0

212.34

2597

16,097

3636.27

2754.80

60.76

60.0

147.06

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

141

Table 8 (Continued) Algorithm

BLS

ILS

MLS

Best-found sequence f(s*)

Fitness values

(1, −z), (29, −z), (28, −z), (41, −z), (40, −z), (19, −x), (23, −x), (52, +x), (50, −x), (48, −x), (56, −x), (62, −x), (55, −x), (54, +x), (60, +x), (17, +x), (21, +x), (16, +x), (53, +x), (47, +x), (39, +x), (58, +x), (59, +x), (20, +x), (57, −x), (63, −x), (18, −x), (22, −x), (61, −x), (9, −x), (7, −x), (8, +x), (6, +x), (2, −y), (3, −y), (5, +y), (4, +y), (46, +x), (51, −x), (49, −x), (38, +x), (10, −z), (42, −y), (43, −y), (30, −z), (31, −z), (35, +y), (13, −z), (14, −z), (11, −z), (12, +y), (15, +y), (24, −z), (32, −z), (25, −z), (26, −z), (36, −z), (37, −z), (27, −z), (45, −z), (44, −z), (34, −z), (33, −z) (40, −x), (1, −x), (28, −z), (41, −z), (20, +x), (62, −x), (29, −z), (16, +x), (55, −x), (17, +x), (22, −x), (19, −x), (57, −x), (54, +x), (52, +x), (60, +x), (50, −x), (61, −x), (18, −x), (53, +x), (63, −x), (56, −x), (59, +x), (58, +x), (21, +x), (47, +x), (39, +x), (48, −x), (23, −x), (8, +x), (6, +x), (2, −y), (51, −x), (9, −x), (7, −x), (3, −y), (5, +y), (46, +x), (38, +x), (4, +y), (49, −x), (10, −z), (42, −y), (43, −y), (30, −z), (35, +y), (31, −z), (14, −z), (13, −z), (11, −z), (24, −z), (32, −z), (25, −z), (12, +y), (15, +y), (27, −z), (36, −z), (37, −z), (44, −z), (45, −z), (26, −z), (34, −z), (33, −z) (57, −z), (48, −z), (1, −z), (41, −z), (29, −z), (40, −z), (28, −z), (55, −x), (23, −x), (52, +x), (47, +x), (17, +x), (18, −x), (61, −x), (60, +x), (59, +x), (54, +x), (16, +x), (20, +x), (21, +x), (53, +x), (50, −x), (63, −x), (39, +x), (58, +x), (56, −x), (62, −x), (22, −x), (19, −x), (8, +x), (6, +x), (5, +y), (49, −x), (9, −x), (3, −y), (38, +x), (2, −y), (4, +y), (7, −x), (30, −z), (46, +x), (42, −y), (10, −z), (13, −z), (51, −x), (31, −z), (11, −z), (35, +y), (43, −y), (14, −z), (24, −z), (12, +y), (32, −z), (25, −z), (36, −z), (15, +y), (37, −z), (26, −z), (27, −z), (33, −z), (44, −z), (34, −z), (45, −z)

Efficiency metrics

fmin

fmax

favg

STD

AvgGap (%)

Success rate (%)

Runtime (s)

2601

6209

3418.50

1211.29

23.76

70.00

218.39

20,162

50,717

32,415.40 11,815.11 1148.67

0.00

235.37

40,060

77,793

58,389.40 12,344.69 2149.21

0.00

231.89

Table 9 The results of sequence planning of the Snap fit electrical connector assembly by SS, GA, GA1, MA, MA1, SA, SA1, BLS, ILS, and MLS methods. Algorithm

Best-found Sequence f(s*)

SS GA GA1 MA MA1 SA SA1 BLS ILS MLS

(8, −y), (1, −y), (3, −y), (5, −y), (4, −y), (6, −y), (2, −z), (7, −z) (8, −y), (1, −y), (3, −y), (4, −y), (6, −y), (5, −y), (2, −z), (7, −z) (8, −y), (1, −y), (4, −y), (3, −y), (6, −y), (5, −y), (2, −z), (7, −z) (8, −y), (1, −y), (3, −y), (4, −y), (5, −y), (6, −y), (2, −z), (7, −z) (8, −y), (1, −y), (4, −y), (3, −y), (5, −y), (6, −y), (2, −z), (7, −z) (8, −y), (1, −y), (4, −y), (6, −y), (3, −y), (5, −y), (2, −z), (7, −z) (8, −y), (1, −y), (4, −y), (6, −y), (3, −y), (5, −y), (2, −z), (7, −z) (8, −y), (1, −y), (3, −y), (4, −y), (6, −y), (5, −y), (2, −z), (7, −z) (8, −y), (1, −y), (3, −y), (4, −y), (5, −y), (6, −y), (2, −z), (7, −z) (8, +z), (1, −y), (4, −y), (3, −y), (6, −y), (5, −y), (2, −z), (7, −z)

Fitness values

Efficiency metrics

fmin

fmax

favg

STD

1018 1018 1018 1018 1018 1018 1018 1018 1018 1019

1018 1020 1921 1021 1020 1020 1021 1019 4250 3073

1018.0 1019.4 1351.5 1019.3 1019.0 1018.6 1018.7 1018.1 2348.2 1757.0

0.00 0.00 0.69 0.14 437.13 32.76 0.95 0.13 0.67 0.10 0.69 0.06 1.06 0.07 0.32 0.00 1330.98 130.67 537.59 72.59

120

0.12

60

100

0.10

50

80

0.08

40

60

0.06

30

40

0.04

20

20

0.02

10

0

0.00 SS MA1 SA

SA1 BLS

AvgGap (%)

Success rate (%)

Runtime (s)

100.0 100.0 60.0 100.0 100.0 100.0 100.0 100.0 40.0 10.0

1.28 5.14 2.52 2.00 1.03 1.28 1.94 1.80 2.28 1.36

0 SS MA1 SA

(a)

(b)

SA1 BLS

SS

MA1

SA

SA1 BLS

(c)

Fig. 12. Average fitness gap with best-known solutions (in percent) of four comparable methods for (a) Tank mockup, (b) Electric connector, and (c) Compact box assemblies.

generating algorithms and computes the relative rank of each value (in case of a tie their average value is set as their common rank) and assigns to its respective algorithm. The sum of the ranks for each algorithm is calculated and through a formula representing the variance of the ranks among algorithms (with

an adjustment for the number of ties), an H statistic is calculated, which approximately has a chi-square distribution. The p-value is then determined as the probability Pr(21 ≥ H), and if it is smaller than the confidence level ˛, then the Null hypothesis is rejected. Numerical results of performing the Kruskal–Wallis test for the

142

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Table 10 The results of sequence planning of the Compact box with snap fit connectors assembly by SS, GA, GA1, MA, MA1, SA, SA1, BLS, ILS, and MLS methods. Algorithm

Best-found Sequence f(s*)

Fitness values fmin

(2, +y), (3, +y), (6, +y), (4, −y), (5, −y), (1, +z), (7, −x), (8, −x), (9, +x), (10, +x), (11, −z), (12, +z), (13, +z), (14, −z), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (1, −y), (2, −y), (4, −y), (5, −y), (6, −z), (10, +x), (9, +x), (11, −z), (3, −z), (8, −x), (7, −x), (15, −x), (12, +z), (13, +z), (14, −z), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (1, −y), (2, −y), (4, −y), (5, −y), (6, −z), (10, +x), (9, +x), (11, −z), (3, −z), (8, −x), (7, −x), (15, −x), (12, +z), (13, +z), (14, −z), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (2, +y), (3, +y), (4, −y), (5, −y), (6, +z), (1, +z), (8, −x), (15, −x), (16, +x), (10, +x), (9, +x), (11, −z), (12, +z), (13, +z), (7, −x), (14, −y), (17, +x), (18, +x), (19, −x), (20, −x) (2, +y), (6, +y), (3, +y), (5, −y), (4, −y), (1, +z), (9, +x), (10, +x), (8, −x), (7, −x), (11, −z), (12, +z), (13, +z), (14, −y), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (2, −y), (5, −y), (4, −y), (3, +y), (6, +y), (1, +z), (11, −z), (8, −x), (7, −x), (9, +x), (10, +x), (12, +z), (13, +z), (14, −y), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (2, +y), (3, +y), (6, +y), (5, −y), (4, −y), (1, +z), (10, +x), (9, +x), (8, −x), (7, −x), (11, −z), (12, +z), (13, +z), (14, −z), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (2, +y), (6, +y), (3, +y), (4, −y), (5, −y), (1, +z), (8, −x), (7, −x), (9, +x), (10, +x), (11, −z), (12, +z), (13, +z), (14, −y), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (4, +y), (2, +y), (5, −y), (6, +y), (3, −y), (1, +z), (11, −z), (10, +x), (12, +z), (8, −x), (7, −x), (9, +x), (13, +z), (14, −y), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (4, −z), (2, −z), (13, +y), (3, +y), (6, +y), (12, +y), (5, −y), (1, +z), (9, +x), (10, +x), (11, −z), (7, −x), (8, −x), (14, −y), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x)

SS

GA

GA1

MA

MA1

SA

SA1

BLS

ILS

MLS

Efficiency metrics

fmax

favg

STD

AvgGap (%)

Runtime (s)

3403.0

3403.0

3403.0

0.00

0.00

100.0

21.93

7003.0

13,284.0

10,242.5

2190.61

200.98

0.0

118.52

7652.6

11,906.0

10,179.2

1386.93

199.12

0.0

132.66

5978.8

11,906.0

8582.3

1775.98

152.19

0.0

136.88

3403.0

5732.4

3806.3

857.73

11.85

80.0

118.84

3403.0

3404.0

3403.1

0.31

0.00

100.0

49.12

3403.0

7402.2

5049.8

1388.60

48.39

20.0

20.14

3403.0

3406.0

3404.4

0.96

0.00

100.0

85.20

4181.8

6631.4

5281.3

691.99

55.19

0.0

20.52

8678.8

13,232.0

11,024.5

1489.14

223.96

0.0

18.73

120

100

100

100

80

80

60

60

40

40

20

20

80

Success rate (%)

60 40 20 0

0

0 SS MA1 SA

SS MA1 SA SA1 BLS

(a)

SA1 BLS

SS MA1 SA

(b)

SA1 BLS

(c)

Fig. 13. Success rates (percentage of reaching feasible solutions) of four comparable methods in 30 runs for (a) Tank mockup, (b) Electric connector, and (c) Compact box assemblies.

Table 11 Ranks and p-values obtained by the Friedman test for all algorithms and benchmark ASP problems. Best values in each row are indicated in bold. Criterion

favg AvgGap (%) Runtime (s) Success rate (%)

Algorithm

p-Value

SS

GA

GA1

MA

MA1

SA

SA1

BLS

ILS

MLS

6.067 5.467 8.400 22.450

23.467 23.633 25.650 12.450

22.217 22.717 24.733 10.450

20.967 21.200 21.133 12.450

12.550 12.583 13.783 19.950

8.583 8.217 11.400 22.450

11.083 10.717 9.667 16.450

9.300 9.033 14.467 20.950

17.650 17.850 15.500 9.450

23.117 23.583 10.267 7.950

0.000 0.000 0.000 0.000

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

143

Table 12 Ratios of mean ranks and p-values obtained by the Kruskal–Wallis test for comparing the SS algorithm with the other nine algorithms in terms of four performance criteria. Entries in boldface indicate outperformance of the SS method. Problem Criterion

1

GA

GA1

MA

MA1

SA

SA1

BLS

ILS

MLS

Rank p

Rank p

Rank p

Rank p

Rank p

Rank p

Rank p

Rank p

Rank p

favg

15.9 45.1

0.0001

15.9 45.1

0.0001

15.9 45.1

0.0001

22.9 38.1

0.0423

23.2 37.8

0.0558

15.9 45.1

0.0001

17.4 43.6

0.0006

15.9 45.1

0.0001

15.9 45.1

0.0001

AvgGap (%)

15.9 45.1

0.0001

15.9 45.1

0.0001

15.9 45.1

0.0001

22.9 38.1

0.0423

23.2 37.8

0.0558

15.9 45.1

0.0001

17.4 43.6

0.0006

15.9 45.1

0.0001

15.9 45.1

0.0001

Runtime (s)

15.9 45.1

0.0002

15.9 45.1

0.0002

15.9 45.1

0.0002

17.9 43.1

0.0012

15.9 45.1

0.0002

15.9 45.1

0.0002

15.9 45.1

0.0002

15.9 45.1

0.0002

15.9 45.1

0.0002

Success rate (%)

45.1 15.9

0.0000

45.1 15.9

0.0000

45.1 15.9

0.0000

34.8 26.2

0.0671

30.5 30.5

1.0000

36.3 24.7

0.0293

34.8 26.2

0.0671

45.1 15.9

0.0000

45.1 15.9

0.0000

favg

17.4 43.6

0.0002

18.9 42.1

0.0006

18.9 42.1

0.0006

18.9 42.1

0.0005

23.2 37.8

0.0124

24.7 36.3

0.0304

29.0 32.0

0.3173

21.8 39.2

0.0051

15.9 45.1

0.0000

(%)

17.4 43.6

0.0002

18.9 42.1

0.0006

18.9 42.1

0.0006

18.9 42.1

0.0005

23.2 37.8

0.0124

24.7 36.3

0.0304

29.0 32.0

0.3173

21.8 39.2

0.0051

15.9 45.1

0.0000

(s)

15.9 45.1

0.0002

15.9 45.1

0.0002

16.2 44.8

0.0002

36.3 24.7

0.1304

30.5 30.5

1.0000

16.5 44.5

0.0002

18.6 42.4

0.0019

15.9 45.1

0.0002

22.9 38.1

0.0534

Success rate (%)

30.5 30.5

1.0000

36.3 24.7

0.0293

30.5 30.5

1.0000

30.5 30.5

1.0000

30.5 30.5

1.0000

30.5 30.5

1.0000

30.5 30.5

1.0000

39.2 21.8

0.0043

43.6 17.4

0.0000

favg

15.9 45.1

0.0000

15.9 45.1

0.0000

15.9 45.1

0.0000

15.9 45.1

0.0000

29.0 32.0

0.3173

18.9 42.1

0.0006

18.9 42.1

0.0006

15.9 45.1

0.0000

15.9 45.1

0.0000

AvgGap (%)

15.9 45.1

0.0000

15.9 45.1

0.0000

15.9 45.1

0.0000

15.9 45.1

0.0000

29.0 32.0

0.3173

18.9 42.1

0.0006

18.9 42.1

0.0006

15.9 45.1

0.0000

15.9 45.1

0.0000

Runtime (s)

15.9 45.1

0.0002

15.9 45.1

0.0002

15.9 45.1

0.0002

15.9 45.1

0.0002

15.9 45.1

0.0002

45.1 15.9

0.0002

21.8 39.2

0.0233

42.1 18.9

0.0025

45.1 15.9

0.0002

Success rate (%)

45.1 15.9

0.0000

45.1 15.9

0.0000

45.1 15.9

0.0000

33.4 27.6

0.1462

30.5 30.5

1.0000

42.1 18.9

0.0004

30.5 30.5

1.0000

45.1 15.9

0.0000

45.1 15.9

0.0000

OSS USS

11 0

AvgGap 2

3

Runtime

12 0

11 0

9 0

ten algorithms, three benchmark problems, and four performance criteria are reported in Table 12. In each column of the table, since all the first three performance criteria (i.e., favg , AvgGap, and Runtime) are of minimization type, ranks smaller than 1 (shown in boldface) indicate that the SS was better than the compared algorithm and criterion for a specific problem. On the other hand, as the fourth criterion (Success rate) is of maximization type, ranks larger than 1 (shown in boldface) indicate the relative outperformance of the SS method. Also, (1 − p) values denote the confidence level of stating that one of the algorithms (e.g. SS in case of bold rank entries) outperforms the other one. The OSS and USS parameters in the last two rows count the times the SS outperformed and underperformed relative to a competing algorithm in a specific criterion. It is interesting to observe that nearly all underperformance cases of the SS relative to other algorithms occurred in the Runtime criterion. The developed SS method proved to be better than other algorithms in all the performance criteria and all three benchmark problems on average, except for the runtimes of SA1, ILS and MLS in the problem 3 (that are respectively 20.14, 20.52, and 18.73 against the 21.93 of the SS as shown in Table 10) which are considerably small and negligible. Note that the SA1, ILS and MLS methods had respectively achieved 20%, 0% and 0% Success rates in producing feasible solutions for the problem 3, which implies that they mostly converged to infeasible solutions in shorter runtimes. Further analysis of the test results against each compared algorithm is as follows: 1. GA, GA1, and MA: The SS algorithm outperformed these algorithms in all criteria for all problems with at least 97.0% confidence, except for the Success rate criterion in the problem 2 in which it was equally fit with the GA and MA algorithms. 2. MA1: For the problem 1, the SS outperformed the MA1 method in all four criteria with at least 93.2% confidence. For the problem 2, although the SS overtook MA1 in the two first criteria, no conclusion can be made about the superiority of either SS or MA1 in the Runtime criterion due to p-values larger than the acceptance level. Also, these algorithms were equally fit in the Success rate criterion for this simple problem. For the problem 3, while the SS outperformed MA1 in the first three criteria (with 99.9% confidence), no conclusion can be made about the superiority of

6 0

3.

4.

5.

6.

10 1

8 0

11 1

11 1

either SS or MA1 in the Runtime criterion due to p-values larger than the acceptance level. SA: The Success rate of the SS algorithm is the same as that of the SA algorithm for all problems. Also these algorithms were equally fit in the Runtime criterion for the problem 2. The SS algorithm dominated the SA in the other criteria for all problems with at least 94.4% confidence except in the first and the second criteria of the problem 3 (the Compact Box), in which the underperformance of the SS still cannot be proved due to the large p-value (0.3173). SA1: The SS algorithm outperformed the SA1 algorithm in all criteria for all problems with at least 97.9% confidence except for the Success rate criterion in the problem 2 in which they were equally fit, and the Runtime of the problem 3 in which the SS algorithm underperformed the SA1 algorithm. The reason is the fast convergence of SA1 to infeasible solutions 80% of the times, which naturally takes less runtimes. BLS: The SS outperformed the BLS method in all four criteria with at least 93.2% confidence for the problems 1 and 3 (except in the success rate criterion for the problem 3). For the problem 2, while the SS overtook BLS in the runtime criterion, the same claim cannot be made with a sufficiently high confidence level for the two first criteria due to the high p-value of 0.3173. ILS and MLS: The SS algorithm outperformed these algorithms in all criteria for all problems with at least 94.6% confidence, except for the runtime criterion in the problem 3 where it took longer runtimes for reaching feasible solutions while ILS and MLS never generated feasible solutions for that problem.

4.2. Flexible vs. rigid ASP analysis In order to show the advantage of considering the deformation behavior of flexible parts in planning assembly sequences, the assembly sequences generated by the SS algorithm for the puzzlelike assembly (Fig. 1) and the three assemblies in Figs. 10 and 11 were solved are compared for both flexible and rigid cases. In the flexible case some parts (as mentioned in the beginning of Section 4) were assumed to be flexible, while in the rigid case all parts were rigid. The results of computational experiments are reported in Table 13 for four benchmark problems in flexible and rigid modes

144

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

Table 13 The results of sequence planning of flexible and rigid assemblies by the SS method. Best values in each column are indicated in bold. Rigidity

Best-found Sequence f(s*)

Fitness values fmin

Soviet medium tank T-34/76 mockup assembly (1, −z), (28, −z), (29, −z), (41, −z), (40, −z), (17, Flexible +x), (52, +x), (16, +x), (20, +x), (54, +x), (60, +x), (53, +x), (39, +x), (58, +x), (21, +x), (59, +x), (47, +x), (18, −x), (56, −x), (50, −x), (19, −x), (55, −x), (61, −x), (48, −x), (62, −x), (23, −x), (57, −x), (22, −x), (63, −x), (9, −x), (7, −x), (8, +x), (6, +x), (3, −y), (2, −y), (5, +y), (4, +y), (38, +x), (46, +x), (51, −x), (49, −x), (10, −z), (31, −z), (30, −z), (43, −y), (42, −y), (35, +y), (13, −z), (14, −z), (11, −z), (24, −z), (32, −z), (25, −z), (12, +y), (15, +y), (26, −z), (36, −z), (37, −z), (27, −z), (44, −z), (34, −z), (45, −z), (33, −z) (1, +y), (29, +y), (28, −z), (41, −z), (40, −z), (39, Rigid +x), (53, +x), (16, +x), (20, +x), (52, +x), (47, +x), (58, +x), (56, −x), (62, −x), (48, −x), (55, −x), (61, −x), (57, −x), (63, −x), (54, +x), (19, −x), (18, −x), (22, −x), (59, +x), (17, +x), (21, +x), (60, +x), (50, −x), (23, −x), (9, −x), (7, −x), (8, +x), (6, +x), (4, +y), (5, +y), (3, −y), (2, −y), (49, −x), (51, −x), (38, +x), (46, +x), (42, −y), (43, −y), (10, −z), (30, −z), (31, −z), (35, +y), (14, −z), (13, −z), (11, −z), (32, −z), (12, +y), (25, −z), (24, −z), (15, +y), (36, −z), (26, −z), (37, −z), (27, −z), (44, −z), (34, −z), (45, −z), (33, −z) Snap fit electrical connector assembly (8, −y), (1, −y), (3, −y), (5, −y), (4, −y), (6, −y), Flexible (2, −z), (7, −z) (1, −y), (8, −y), (5, −y), (3, −y), (6, −y), (4, −y), Rigid (2, −z), (7, −z) Compact box with snap fit connectors assembly Flexible (2, +y), (3, +y), (6, +y), (4, −y), (5, −y), (1, +z), (7, −x), (8, −x), (9, +x), (10, +x), (11, −z), (12, +z), (13, +z), (14, −z), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) (1, +z), (9, −z), (5, −z), (8, −z), (4, −z), (10, −z), Rigid (2, −z), (6, −z), (7, −z), (3, −z), (11, −z), (12, +z), (13, +z), (14, −y), (15, −x), (16, +x), (17, +x), (18, +x), (19, −x), (20, −x) Puzzle-like assembly Flexible (3, −x), (1, −x), (4, −x), (5, −x), (6, −x), (2, −x) Rigid (3, −x), (4, −x), (5, −x), (6, −x), (2, −x), (1, −x)

Efficiency metrics fmax

favg

2596

2602

9020

33,334

1018

AvgGap (%)

Success rate (%)

Runtime (s)

1.84

0.18

100.0

139.79

24,778.7

8544.19

854.49

0.0

220.82

1018

1018.0

0.00

0.00

100.0

1.28

2701

2701

2701.0

0.00

165.32

0.0

0.46

3403

3403

3403.0

0.00

0.00

100.0

21.93

10,806

13,507

11,526.5

1022.25

238.72

0.0

12.55

441 900

441 901

441.0 900.1

0.00 0.32

0.00 104.10

100.0 0.0

and fmin , fmax , favg , AvgGap (%), Runtime (s), and Success rate (%) performance criteria. The results showed that incorporating flexibility in the ASP model resulted in better results in all performance criteria for all four benchmark problems on average compared to mere rigid modeling, except in runtimes of problems 2 and 3 (with 0.46 and 12.55 seconds against 1.28 and 21.93 seconds, respectively), which are not significant. Also, while the flexible ASP attained 100% success rates, the rigid ASP failed to produce feasible solutions for all benchmark problems. This shows that considering flexibility of parts helps to obtain feasible assembly sequences for parts that are unassemblable when considered rigid. The feasibility of a solution is checked by examining if the required assembly stresses (i.e., entries in the ASM) of each part and all parts assembled prior to it (i.e., those appearing before it in the solution string) does not exceed  m as in (4). In order to statistically approve the significance of the outperformance of flexible over rigid ASP, a “Friedman two-way ANOVA by ranks” and a “Kruskal–Wallis one-way ANOVA by ranks” test were conducted for the four benchmark problems, the results of which are reported in Tables 14 and 15. The p-values in both tables show that the performance of the SS algorithm with flexible parts in favg , AvgGap, and Success Rate criteria was decisively superior to when

2600.60

STD

0.12 0.13

all parts are considered rigid (with a confidence level of 100.0%). However, in the Runtime criterion for the electrical connector and compact box assemblies, the rigid ASP outperformed the flexible ASP since it converged to infeasible local minimum solutions in shorter times while the flexible ASP spent more time on searching deeper for feasible solutions. The worst-case time complexity order of constructing the ASM equals to that of constructing the AIM. Both rigid and flexible ASP methods require formation of the AIM, in which the statuses of all pair-wise interferences between parts are recorded. The time needed for checking all possible collisions between any two parts is proportional to p(n, 2) = n(n − 1), which is the number of permutations of all pairs of parts (note that combinations are not considered since pair-wise blockages are not symmetric). Since in the AIM, blockages must be checked along 2 m directions, n(n − 1)2 m entries must be calculated, which is in the order of O(mn2 ). For constructing the ASM, any pi , pj , and dk for which Ii,k j = 1 in the AIM and one (or both) of the parts are flexible, their corresponding assembly operation must be analyzed via FEM (in constant time) to check if the flexibility of the parts resolves the blockage. In the worst case, when all parts are flexible, all elements of the AIM matrix (except for those on the main diagonal) must be analyzed in the FEM software, and so the order of forming the ASM is

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

145

Table 14 Ranks and p-values obtained by the Friedman test for flexible and rigid benchmark ASP problems. Criterion

favg

AvgGap (%)

Runtime (s)

Success rate (%)

Flexible Rigid

2.600 6.400

2.500 6.500

4.650 4.350

6.500 2.500

p-value

0.000

0.000

0.578

0.000

Table 15 Ratios of mean ranks and p-values obtained by the Kruskal–Wallis test for comparing flexible and rigid ASP in terms of four performance criteria. Entries in boldface indicate outperformance of the Flexible vs. rigid ASP. Problem

Electrical connector

Compact box

Rank

Tank mockup p

Rank

p

Rank

p

Rank

p

Criterion favg

15.9 45.1

0.0001

15.9 45.1

0.0000

15.9 45.1

0.0000

15.9 45.1

0.0000

AvgGap (%)

15.9 45.1

0.0001

15.9 45.1

0.0000

15.9 45.1

0.0000

15.9 45.1

0.0000

Runtime (s)

15.9 45.1

0.0002

45.1 15.9

0.0001

45.1 15.9

0.0002

21.8 39.2

0.0165

Success rate (%)

45.1 15.9

0.0000

45.1 15.9

0.0000

45.1 15.9

0.0000

45.1 15.9

0.0000

O(mn2 ). The time required for finite element analysis of a pair of parts mainly depends on the type of meshing and the number of their elements, which for our benchmark assemblies was in the order of few minutes. However, such extra preprocessing times have not been considered in the comparative statistical tests as we wanted to evaluate the impact of using ASM (that is, considering the parts’ flexibility and hence assemblability of the product thereof) in the assembly sequence planning algorithm instead of the simpler AIM which just implies geometrical blockages of rigid parts. 5. Conclusions Assembly sequence planning (ASP) of products with both rigid and flexible parts is a challenging problem on the way of automatizing assembly processes of many real-world products. Thus far, the ASP for rigid parts have been widely addressed in the literature, but deformability conditions of flexible parts has not been taken into account. This paper lays out a theoretical ground for modeling the deformability of flexible assembly parts by introducing the Assembly stress matrix (ASM), and a method is developed to generate and optimize assembly sequences for products with rigid and flexible parts, to the best of our knowledge, for the first time. This is done by simultaneously minimizing two operational criteria of total maximum applied assembly stress and number of assembly direction changes. At first, a new heuristic procedure is implemented for generating a proper initial population which is improved by a novel Effective Insertion operator that tries to make the solutions in the initial population as feasible as possible. The objective function used for evaluating the quality of an assembly sequence is a weighted sum of the number of direction changes and the required force to perform the assembly operations. The algorithm used for optimizing this objective function is the powerful Scatter Search (SS) method with its parameters tuned by a TOPSIS-Taguchi-based method that considers the following criteria: (1) number of fitness function evaluations, (2) average total maximum applied assembly stress, (3) average number of same-direction assembly operations in one assembly sequence, and (4) the used objective function. Three ASP problems selected from the literature or designed by us were solved by the presented SS and 9 other metaheuristic methods, while the parameters of all competing algorithms were tuned by the same TOPSISTaguchi-based parameter tuning method. Computational results showed that the SS outperformed the other algorithms by producing the best-known or optimal solutions most of the time thanks

Puzzle-like assembly

to the careful selection of the initial population, reference set selecting and updating, subset generation, solution recombination, and improvement method. The outperformance of the SS method is verified by conducting Friedman two-way and Kruskal–Wallis one-way ANOVA tests per each product, which approved the significance of its lead over other algorithms. Also, an in-depth analysis of flexible vs. rigid ASP has been done which shows that considering the deformation behavior of flexible parts can lead to finding feasible assembly sequences for products that are unassemblable with rigid parts. This claim was also verified by conducting statistical tests for four products in both flexible and rigid modes. As a result, based on the computational results, and statistical tests, the developed SS-based ASP solver can help design and manufacturing engineers in finding higher quality assembly sequences, which in turn reduces the product’s overall production costs. In the current work we have assumed that the ASP problem is sequential, linear, and monotone: meaning that the assembly can be accomplished by two hands only, parts are added to the assembly separately (and not as a subassembly), and each part is manipulated just once, directly from its initial position to its final, desired configuration (see [18] for definitions). However, while the majority of relevant works in the literature consider the above assumptions commonly, an open problem, which we believe is very challenging, is to solve non-sequential, nonlinear, and non-monotone ASP problems with flexible parts. Another issue with the present and many other works is that only straight (main Cartesian) assembly directions are considered, while many real-world assemblies require complex motions for assembling parts. So a challenging problem is to consider rotations, curved trajectories, and more complex deformations such as torsion, shear, bend, and fold during the assembly of flexible parts. Another issue that has not been considered in this paper is the parts’ tolerances and variations in size (e.g. deviations from nominal dimensions such as length and diameter), in form (e.g. deviations from the nominal shape such as out-of-round holes), in surface (e.g. small irregularities in surfaces due to the manufacturing process), and in kinematics (which occur when mating parts adjust their positions slightly due to gaps or interferences). However, we believe that the proposed method can still be applied to toleranced parts by setting the required assembly stress entries in the ASM for the toleranced parts equal to the sum of their Yield stresses d (e.g. S = Yi + Yj ), since yield stresses are the maximum stresses i i ,j applicable on parts without inflicting plastic deformation. While through this approach the fitness values of all generated solutions may become a bit larger than their real values, since this condition is

146

S. Ghandi, E. Masehian / Journal of Manufacturing Systems 36 (2015) 128–146

established for all generated solutions, still the best (global optimal) solution remains unchanged, which it is again the solution with the smallest amount of the objective function. Nevertheless, assemblability of toleranced parts is an extensive and challenging research topic since the required forces/stresses for assembling such parts could be significantly different for one part from another, with risks of causing plastic deformation when assembling transition or interference fits. References [1] Acker J, Henrich D. Manipulation of deformable linear objects: from geometric model towards program generation. In: Proceedings of the IEEE international conference on robotics and automation. 2005. p. 1541–7. [2] Auray D., Kiely K.M. Snap fit electrical connector assembly with operating tool for facilitating the connection of a connector assembly to an electrical box; 2007. US Patent No. 07205489. Available http://www.strutpatent.com/patent/07205489/snap-fit-electricalat: connector-assembly-with-operating-tool-for-facilitating-the-connection-ofa-connector-assembly-to-an-electrical-box [3] Barr AH. Global and local deformations of solid primitives. ACM Siggraph Comput Graph – ACM 1984:21–30. [4] Bonneville F, Perrard C, Henrioud J-M. A genetic algorithm to generate and evaluate assembly plans. In: Proceedings of INRIA/IEEE symposium on emerging technologies and factory automation. 1995. p. 231–9. [5] Cao P-B, Xiao R-B. Assembly planning using a novel immune approach. Int J Adv Manuf Technol 2007;31:770–82. [6] Chen CP. Design of a real-time AND/OR assembly scheduler on an optimization neural network. J Intell Manuf 1992;3:251–61. [7] Chen W-C, Tai P-H, Deng W-J, Hsieh L-F. A three-stage integrated approach for assembly sequence planning using neural networks. Expert Syst Appl 2008;34:1777–86. [8] Ching-Lai H, Yoon K. Multiple attribute decision making: methods and applications. Springer-Verlag; 1981. [9] Derrac J, García S, Molina D, Herrera F. A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm Evol Comput 2011;1: 3–18. [10] Friedman M. The use of ranks to avoid the assumption of normality implicit in the analysis of variance. J Am Stat Assoc 1937;32:675–701. [11] Gao L, Qian W, Li X, Wang J. Application of memetic algorithm in assembly sequence planning. Int J Adv Manuf Technol 2010;49:1175–84. [12] Gibson SF, Mirtich B. A survey of deformable modeling in computer graphics. Technical Report TR97-19. Cambridge, MA: Mitsubishi Electric Research Laboratories; 1997. November 1997. Available at: http://www.merl.com/ publications/TR97-19 [13] Glover F. Heuristics for integer programming using surrogate constraints. Decis Sci 1977;8:156–66. [14] Glover F. Tabu search and adaptive memory programming – advances, applications and challenges. In: Interfaces in computer science and operations research. Springer; 1997. [15] Guo J, Wang P, Cui N. Adaptive ant colony algorithm for on-orbit assembly planning. In: 2nd IEEE conference on industrial electronics and applications. 2007. p. 1590–3. [16] Hong D, Cho H. A neural-network-based computational scheme for generating optimized robotic assembly sequences. Eng Appl Artif Intell 1995;8: 129–45. [17] Hong D, Cho H. Generation of robotic assembly sequences using a simulated annealing. In: Proceedings of IEEE/RSJ international conference on intelligent robots and systems. 1999. p. 1247–52. [18] Jiménez P. Survey on assembly sequencing: a combinatorial and geometrical perspective. J Intell Manuf 2013;24:235–50.

[19] Kruskal WH, Wallis WA. Use of ranks in one-criterion variance analysis. J Am Stat Assoc 1952;47:583–621. [20] Laguna M, Marti R, Martí RC. Scatter search: methodology and implementations in C. Springer; 2003. [21] Luo Q, Xiao J. Haptic rendering involving an elastic tube for assembly simulations. In: The 6th IEEE international symposium on assembly and task planning: from nano to macro assembly and manufacturing. 2005. p. 53–9. [22] Lv H, Lu C. An assembly sequence planning approach with a discrete particle swarm optimization algorithm. Int J Adv Manuf Technol 2010;50:761–70. [23] Martí R, Laguna M, Campos V. Scatter search vs. genetic algorithms. In: Metaheuristic optimization via memory and evolution. Springer; 2005. [24] Masehian E, Akbaripour H, Mohabbati-Kalejahi N. Landscape analysis and efficient metaheuristics for solving the n-queens problem. Comput Optim Appl 2013;56:735–64. [25] Motavalli S, Islam A-U. Multi-criteria assembly sequencing. Comput Ind Eng 1997;32:743–51. [26] Ranjit R. A primer on the Taguchi method. van Nostrand Reinhold; 1990. [27] Rashid MFF, Hutabarat W, Tiwari A. A review on assembly sequence planning and assembly line balancing optimisation using soft computing approaches. Int J Adv Manuf Technol 2012;59:335–49. [28] Remde A, Henrich D, Wom H. Manipulating deformable linear objects-Contact state transitions and transition conditions. In: Proceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems. 1999. p. 1450–5. [29] Sacerdoti ED. Planning in a hierarchy of abstraction spaces. Artif Intell 1974;5:115–35. [30] Sederberg TW, Parry SR. Free-form deformation of solid geometric models. ACM Siggraph Comput Graph – ACM 1986:151–60. [31] Sinanoglu C, Börklü HR. An assembly sequence-planning system for mechanical parts using neural network. Assem Autom 2005;25:38–52. [32] Su Q. Computer aided geometric feasible assembly sequence planning and optimizing. Int J Adv Manuf Technol 2007;33:48–57. [33] Talbi E-G. Metaheuristics: from design to implementation. John Wiley & Sons; 2009. [34] Tiwari M, Kumar A, Mileham A. Determination of an optimal assembly sequence using the psychoclonal algorithm. Proc Inst Mech Eng Part B: J Eng Manuf 2005;219:137–49. [35] Tong L-I, Wang C-H, Chen H-C. Optimization of multiple responses using principal component analysis and technique for order preference by similarity to ideal solution. Int J Adv Manuf Technol 2005;27:407–14. [36] Tseng Y-J, Yu F-Y, Huang F-Y. A green assembly sequence planning model with a closed-loop assembly and disassembly sequence planning using a particle swarm optimization method. Int J Adv Manuf Technol 2011;57:1183–97. [37] Wang J, Liu J, Zhong Y. A novel ant colony algorithm for assembly sequence planning. Int J Adv Manuf Technol 2005;25:1137–43. [38] Wang L, Hou Y, Li X, Sun S. An enhanced harmony search algorithm for assembly sequence planning. Int J Model Identif Control 2013;18:18–25. [39] Wang Y, Liu J. Chaotic particle swarm optimization for assembly sequence planning. Robot Comput-Integr Manuf 2010;26:212–22. [40] Wolter J, Kroll E. Toward assembly sequence planning with flexible parts. In: Proceedings of IEEE International Conference on Robotics and Automation. 1996. p. 1517–24. [41] Zeng B, Li MF, Zhang Y. Research on assembly sequence planning based on firefly algorithm. J Mech Eng 2013;11:025. [42] Zhang H, Liu H, Li L. Research on a kind of assembly sequence planning based on immune algorithm and particle swarm optimization algorithm. Int J Adv Manuf Technol 2013:1–14. [43] Zhou W, Yan J, Li Y, Xia C, Zheng J. Imperialist competitive algorithm for assembly sequence planning. Int J Adv Manuf Technol 2013;67:2207–16. [44] Zhou W, Zheng J-R, Yan J-J, Wang J-F. A novel hybrid algorithm for assembly sequence planning combining bacterial chemotaxis with genetic algorithm. Int J Adv Manuf Technol 2011;52:715–24. [45] Zvezda C. Миhиɑtюры-OOO Звeздɑ; 2014. Available at http://www.zvezda. org.ru/?lng=1&nav=2&p=59&set=6101. Also at http://www.hobbyandleisure. co.uk/hlstore/catalog/zvezda-5001-soviet-medium-tank-3476-mod-1943172-t34-p-6660.html