Optimization of Bearing Locations for Maximizing First Mode Natural Frequency of Motorized Spindle-Bearing Systems Using a Genetic Algorithm

This paper has developed a genetic algorithm (GA) optimization approach to search for the optimal locations to install bearings on the motorized spindle shaft to maximize its first-mode natural frequency (FMNF). First, a finite element method (FEM) dynamic model of the spindle-bearing system is formulated, and by solving the eigenvalue problem derived from the equations of motion, the natural frequencies of the spindle system can be acquired. Next, the mathematical model is built, which includes the objective function to maximize FMNF and the constraints to limit the locations of the bearings with respect to the geometrical boundaries of the segments they located and the spacings between adjacent bearings. Then, the Sequential Decoding Process (SDP) GA is designed to accommodate the dependent characteristics of the constraints in the mathematical model. To verify the proposed SDP-GA optimization approach, a four-bearing installation optimazation problem of an illustrative spindle system is investigated. The results show that the SDP-GA provides well convergence for the optimization searching process. By applying design of experiments and analysis of variance, the optimal values of GA parameters are determined under a certain number restriction in executing the eigenvalue calculation subroutine. A linear regression equation is derived also to estimate necessary calculation efforts with respect to the specific quality of the optimization solution. From the results of this illustrative example, we can conclude that the proposed SDP-GA optimization approach is effective and efficient.


Introduction
Spindle-bearing systems have been widely applied in mechanisms, such as machine tools, that need relative rotary motion to accomplish desired machining functions.Recently, to alleviate high-speed rotation, motorized spindles have been successfully developed and used in high-speed machine tools.This type of spindle is equipped with a built-in motor as an integrated part of the spindle shaft, eliminating the need for conventional powertransmission devices.To achieve high cutting performance, the motorized spindle-bearing system must be designed deliberately since it is the main moving component in the machine tool.
The design process for a motorized spindle system usually begins with identification of the bearing types, configurations, preloading methods, and spindle motor systems that would be appropriate for the desired machining application [1].Once the bearings and spindle system concept has been clarified, detailed product specifications are further used to define other design values of the spindle system.One example is the locations of bearings assembled on the spindle shaft.However, in the initial stage of the spindle-bearing system design, the engineers must ascertain the optimal values of design variables that maximize the first-mode natural frequency (FMNF) of the spindle-bearing system to avoid forced resonance from easily happening in machining operations [2].
How the spindle shaft and bearing related design variables affect a system's static and dynamic characteristics has been intensively studied [3]- [7].Recently, Lin and Tu [2] comprehensively explored the effects of design variables on the natural frequencies of a high-speed spindle system.The eight design variables they considered included material of spindle shaft, diameters and total length of the spindle shaft, bearing initial preload, spacings between the bearings of the front or rear bearing set, spacing between the middle line of the front and rear bearing sets, and spacing of the middle line of the front and rear bearing sets to the end of the cutting tool.Their results showed that the first two most important design factors among the eight were related to the locations of the bearings, which implied that the positions bearings installed on the spindle shaft must be determined carefully.There have been several papers dedicated to the optimization of locations or spans of bearings based on the stability of rotating machinery or machining from the practical viewpoints of machining experts and scholars [6]- [8].However, there still lacks a research concentrated directly on the dynamics of spindle bearing systems.
To estimate natural frequencies of spindle systems analytically in the early detailed design stage, the finite element method (FEM) has been frequently adopted in modeling rotor dynamics, due to its flexibility in treating ever-complex spindle system designs, such as motorized spindles.Basically, the FEM model for the spindle systems of machine tools is similar to those developed in rotor dynamic literature [9]- [12].However, the spindle shafts used in machine tools usually have smaller shaft diameters and bearings, and possess no disk-like components.As the FEM model is built, the system natural frequencies can then be solved as an eigenvalue problem of the dynamic model, which is equal to finding the roots of a high order polynomial equation [13]; therefore, the FMNF maximization problem generally appears nonlinear.Because derivatives of the objective function FMNF are difficult or even impossible to deduce, genetic algorithm (GA), an optimization approach without derivatives, is utilized to solve the optimization problem in this paper.
Genetic algorithm, first introduced by Holland in 1975 [14], is a modern metaheuristic method which can solve nonlinear optimization problems effectively.GA is a computational technique which simulates the evolution process based on the principles of natural genetics and natural selection [15].In the academic field of engineering optimization, GA has been successfully applied in the researches such as identification of a spatial slider-crank mechanism [16], optimal designs of pressure swing adsorption [17], comparisons for damage detection on structures [18], and rolling element bearing design [19].
The main purpose of this paper is to demonstrate the application of GA on determining the optimal locations of bearings of motorized spindle shafts, which are the decisions expected to be made by design engineers when they start to plan assembling the bearings and spindle shaft after the spindle and bearing specifications have been specified during the concept design stage.In the following sections, the FEM model of a typical motorized spindle system is introduced first.Next, based on the FEM model, an optimization mathematical model constituted by an objective function and constraints represented with the design variables, i.e. the locations of bearings, is developed, and an appropriate GA-based optimal solutions searching approach is described.Then, to manifest the proposed GA, a spindle-bearing system design problem is illustrated.The mathematical model of the bearing location problem for the illustrative design is built first.To find the optimal solutions, a customized computer code is written in MATLAB by following the procedures of the proposed GA, where the GA parameters are determined by statistical analyses.The computer program is also used to investigate the influence of the parameters on the results.Finally, the primary conclusions and discussions of this paper are described.

FEM Model of Spindle-Bearing Systems
Without loss of generality, the proposed dynamic FEM model is developed based on the spindle-bearing system as shown in Figure 1, where two sets of angular-contact ball bearings (Set 1 and Set 2) are utilized.Figure 2 shows the essential mechanical models and the major independent design variables required to describe the dynamic properties of spindle-bearing systems [2], in which the integrated spindle-bearing system FEM model is a combination of the distributed spindle shaft FEM and the discrete bearing stiffness models.Those models can be represented as functions of spindle specifications, bearing specifications, initial preload, and bearing location.In this research, however, only the bearing location (highlighted in Figure 2) is retained as the design variables while the other factors are treated as design parameters.
To construct the spindle shaft FEM model, the spindle shaft is first discretized into a finite number of beam elements, and each node of the elements is assigned with four degrees of freedom, where, in sequence, the first two are assigned to the lateral directions and the remaining two are assigned to the angular directions.Associated with several specific shape functions, the kinetic and potential energy of each element can be obtained by integrating those of the cross-sections along the axis of the spindle shaft, and expressed as functions of the physical and geometrical properties of the element.Summing up the kinetic and potential energy of all elements, we can get the total energy of the spindle shaft, and by utilizing Lagrange's equation, the equations of motion (EOM) of the spindle shaft can be finally deduced.After combining the bearing stiffness with the EOM of the spindle shaft, the simple free vibration EOM for the motorized spindle-bearing system can be presented as [20] 0 where q in ( ) is the global node displacement vector of the spindle shaft and p is the element number; M is the mass matrix of size ( ) ( )  where b n is the total number of bearings in the system and i B , with the same size as s K , is the stiffness matrix contributed by bearing i .Usually the points of application of the bearings are arranged on the joint stations of spindle shaft FEM model [13].If we neglect the axial stiffness contributed by the bearing as in the analysis of the spindle shaft, the entries of θ [21]: where a c is an empirical data decided by the experimental results.In order to find the natural frequencies of the system, we relate the vibration problem to the generalized eigenvalue problem by substituting e j t ω = q h into Equation (1), which results in [13] ( ) where ω is the system natural frequency, λ is the eigenvalue, and h is the eigenvector.
To find a nonzero solution h of Equation ( 6), we must have [22] ( ) which is the characteristic equation of

M K are
( ) , and therefore, the natural frequencies are ( ) , , , for each eigenvalue, there is a corresponding eigenvector.

Mathematical Model of the Bearing Location Optimization
In this research, we define the locations of bearings as the design variables.Assume that there are s n sets of bearings to be installed on the spindle system, and for set i , i n identical bearings are used and, therefore, Under this condition, the locations of bearings of set i can be represented as , , , where the real number , i j x  , named local design variable, indicates the coordinate of the middle line of bearing j of the bearing set i .The lower the bearing sets numbered, the closer their locations are to the original point.Since the bearings are identical within a set, index j can also be used to represent the sequence of the bearings such that the bearings with lower j 's are closer to the original point.Therefore, for bearing set i , ,1 i x  and are the coordinates of the two end bearings respectively.The adjacent bearings of bearing j are bearings 1 j − and 1 j + for 2,3, , 1 . The locations of all bearings x can be formed by assembling all i Besides the local variables, we define the global design variables as , which are utilized in calculating natural frequencies after being converted into the corresponding numbers of joint stations in FEM, and for a local variable owning an index ( ) , i j , the index k for its representing global variable can be evaluated by ( ) with 0 0 n = .Since the FMNF is the smallest among the natural frequencies of spindle-bearing systems, maximizing it draws the following objective function 1 Max ω (12) where 1 min ω = ω .As discussed in Section 2, ω is the square root of eigenvalue λ and depends on x since the matrix K is a function of x .
There are two different kinds of location-related constraints for the design variables considered in this research.First, the locations of bearings are ordinarily constrained by, for example, the geometrical boundaries of the segments where they are located.This kind of constraint is set to avoid interference with other parts, and or by the requirements originated from the results of other analyses such as statics.Here we name those constraints as Type I constraints and represent them with , , , , , where , i j L and , i j U are the upper and lower limits of , i j x  respectively.The second type of constraints, Type II constraints, reflects the fact that the least amount of spacing between two bearings must be greater than a specified value such as the width of the bearing.Apparently, only the bearings within the same group would be involved in the same constraint equations of Type II constraints.If the smallest spacing of bearings j and 1 j + in the bearing set i is recognized as , , 1 i j j C + , the constraints can be written as , , 1 Ideally, the optimal solutions are obtained by identifying the ones possessing the maximum FMNF from the feasible solutions.However, since the objective function is complicated and nonlinear, it may be difficult or even impossible to derive an analytical or gradient-based method to solve the optimization problem.In this paper, a GA-based approach is constructed instead to search for the optimal solutions, which simultaneously satisfy all the constraints presented in Equation ( 13) and Equation ( 14).

Formulation of the Genetic Algorithm
As a randomized population-based search technique, GA utilizes genetic operators, i.e. crossover and mutation, inspired by natural evolution process to manipulate individuals in a population over generations to improve their fitness gradually based on the "survival of the fittest" strategies.The individuals in GA are likened to chromosomes which are most commonly represented as strings with an equal length of binary numbers formed by the coded design variables, and to each chromosome, there corresponds a value of the objective function, referred to as the fitness of the chromosome [23].
To encode the real variables x as binary numbers for forming the chromosome, we adopt the simple binary representation scheme with m bits [24].The value of bit number m is generally determined based on the required precision of the design variables.For each component i x , by translating and scaling, we map its fea- sible region , ssed as m -bit binary strings i y , which defines the corresponding encoding variable of i x as ( ) ( ) where  , and by collecting all i y 's, we can form , the design variables in binary format.As a result, the chromosome S can be constituted by attaching the m -bit binary strings of all variables, end for end, in sequence as S y y y a a a a a a a a a where c b n n m = × is the length of chromosome S .S can be converted back to y with a similar reverse transformation.
When we perform the decoding tasks to convert the binary codes back to real numbers, the feasible regions of variables must be identified deliberately according to the constraints applied to them.For a simple problem with fixed upper and lower limits, we can conduct decoding on all variables simultaneously since they are independent to each other.However, for mutually dependent variables such as those considered in this paper, their values in new generations, yielded by genetic operations of GA, may violate the dependency constraints (such as Type II constraints) if we try to decode the variables parallel in time or fail to handle the dependencies of the variables properly.There are several methods to ensure variable dependent constraints satisfied during GA optimization.The most effective approach is to restrict the search to valid regions of the search space, i.e. the chromosome is decoded in such a way that invalid solutions are prevented.We can thereby avoid wasting effort by evaluating infeasible solutions [25].
In this research, to ascertain that the design variables satisfy all constraints after decoding from binary codes to real values, we propose the Sequential Decoding Process (SDP), in which the valid regions of variables are decided one by one, sequentially.Since the locations of bearings of two different sets are independent, we can simultaneously decode the variables within different bearing sets at the same time, and for each set of bearings, we decode the variables, from the first to the last, in sequence.The first step of SDP is converting the global binary variable y to the local binary variable , i j y  by using Equation (11).Then, for each bearing set i , 1, 2, , s i n =  , the binary , i j y  , with respect to global k y , is decoded into real , i j x  by using the standard binary decoding formula for j from 1 to i n , in sequence: where if 1 j = , the valid region , , , , , ,  11).Now that the above representation scheme has been defined, the procedure of applying the proposed SDP and GA, which is called SDP-GA, to search for the optimal solutions of the mathematical model is shown in Figure 3 and summarized as [26]: 1) Utilize a random number generator to create an initial population which consists of z individuals, each with an c n binary string.3) Create a mating pool by using the proportional selection procedure, such as a simulated roulette wheel, to choose the members in the initial population to form a mating pool with a size z , in which the chance of an individual being selected for the mating pool is proportional to its fitness value.4) Perform crossover operation on the mating parent pool to produce offspring.The crossover operation adopted in this research involves randomly choosing two chromosomes from the mating pool and an integer k in the range of 1 to 1 c n − .The first k positions of the parents are exchanged to produce two children.5) Perform mutation operation by visiting every bit of all individuals of the new population and switching the bit from 0 to 1 or 1 to 0 with a mutation probability m p . 6) Evaluate the population as what we did in step 2. The highest fitness value and the corresponding variable values are stored.A generation is now completed.7) If the preset number of generations g n is complete, the process is stopped; otherwise, go to step 3. Note that the size of population z , the number of generations g n , and the mutation probability m p are ge- nerally decided with a few experiments to ensure a good convergence.However, to further determine the optimal values of the parameters in this research, we conduct statistical analysis such as design of experiments (DOE), analysis of variance (ANOVA), and linear regression [27].

An Illustrative Example
To demonstrate the validity of the proposed SDP-GA, a real bearing location optimization problem of a high-speed spindle-bearing system design is exemplified in this section.It is assumed that the concept design stage has already been completed, and the results are summarized as specifications in Table 1 and Table 2 for the spindle shaft and bearings, respectively.The spindle shaft in this illustrative case is made of steel and constitutes 14 cylinders with different diameters.Table 1 shows the material properties of the spindle shaft, and the length and diameter of each segment as well as the number of elements that each segment is divided into in the FEM model.There are two different kinds of bearings, types ACB1 and ACB2, used in the example, and their specifications are indicated in Table 2.The table also includes the widths of the bearings and the data required in Equation ( 5) for calculating their stiffnesses with the assumption that preloads applied on bearings are not affected by their locations.The results are also provided in the last column of the table.
If in the system level design stage, four bearings, two of type ACB1 on segment 7 and two of type ACB2 on 12, are planned to be installed in the system, i.e.    , , , where the set formed by the two bearings on segment 7 is numbered one and the set on segment 12 is numbered two.Our objective is to maximize 1 ω of the spindle system.The constraints of the design variables can be derived from the data provided in Table 1 and Table 2. Assumed that geometry limits are the only sources of Type I constraints.These constraints for the design variables are illustrated as Figure 4.As revealed in Figure 4(a), for example, the lower limit 1,1 L of 1,1 x  can be obtained by adding the left end coordinate of segment 7 and half of the width of the type ACB1 bearing, i.e. 1,1 155 0.5 18 164 U is equal to subtracting the right end coordinate of segment 7 off 1.5 times of the width of type ACB1 bearing, i.e. x  and 1,2 x  is the width of the bearing ACB1, which is equal to 18 mm, and between 2,1 x  and 2,2 x  is 16 mm, the width of bearing ACB2.Table 3 summarizes all constraints for the four local design variables and also their valid regions utilized in SDP-GA.
Since there are only four variables in this case, we can reach the global optimal solution as [ ] 189, 296,558,574 * = x mm and the corresponding 1 ω * being 794.622Hz by an exhaustive search with the precision as 0.5 mm.However, even with only four variables, to obtain the results, the FMNF-related eigen-  x  times, which takes about 240 days if one execution takes 0.05 seconds.To a design project, 240 days may be too long to wait for a decision and the schedule is likely to be delayed.To avoid such a tedious, exhaustive search, the alternative GA approach provides a more efficient approach in obtaining an acceptable, near optimal solution.
To code a computer program for searching for the optimal solutions by following the SDP-GA procedures developed in Section 4, we need to decide bit numbers of the design variables first.If the precision is required to be at least 0.5 mm and since the largest interval for the design variables is 114 mm, the binary number m , in order to guarantee the required precision, must satisfy ( ) 114 2 1 0.5 m − ≤ , which leads to 8 m = , and such that the length of chromosome 32 c n = .Following the procedures indicated in Figure 3, a computer program which constitutes a main script M-file to execute the SDP-GA and a supporting function M-file to calculate the FMNF, is coded in MATLAB.However, we still need to determine the parameters z , g n , and m p before running SDP-GA in searching for the optimal solutions.
In this paper, we use the statistics techniques DOE to identify the optimal values of the parameters under the practical assumption that the allowable time to decide the locations of bearings is one minute, i.e. the number of calculations for finding the eigenvalues must be 1200 at most if each eigenvalue evaluation takes 0.05 seconds.As the total number of eigenvalue calculations ( ) in the SDP-GA optimization process, to satisfy the required calculation number 1200, possible values of the parameters z and g n can be expressed as ( ) , g z n combinations that include (16,74), (20,59), (30,39), and (40,29), where the smallest z is set as 16 since it is the smallest number of chromosomes to prevent a population from ending up with identical chromosomes.Because z and g n are dependent, we use the levels of z to represent the combinations.The third parameter m p is typically very small [24], therefore, the possible values considered here are 0.005, 0.01, and 0.02.
To examine convergence, Figure 6 shows the evolution of the optimization process in a certain test run for all four candidates ( ) , g z n with 0.01 m p = respectively, where the maximum and average values of FMNF for each of the generations are presented.From the figures, we can find that, as generations go by, not only does the best chromosome of the population improve, but the rest in the entire population also improve.This implies that the proposed SDP-GA approach not only provides a single solution, but also a family of good-quality solutions.
After preliminarily confirming the convergence of the proposed SDP-GA algorithm, to decide the levels for parameters z and m p , a two-factor experiment is performed with FMNF as the response variable and with z and m p as the main factors.Since z and m p are with 4 and 3 levels respectively, there are 12 different combinations or cells of the two factors.In a complete balanced experimental design, we obtain 30 observations (replication measurements) from each of the 12 cells.The sums, sums of squares, averages, and standard deviations of the 30 observations for each cell are calculated and shown in = , there is no indication of interaction between these two factors.The last column of Table 5 shows that the p -value for the test statistic for the main effect z is considerably less than 0.05, while for the main effect m p and the interaction are greater than 0.05.A graph of the FMNF averages compared to levels of z for each m p is shown in Figure 7. Notice that of the 12 experimental configurations considered, the optimum configuration employs 16 z = at 0.01    figure, where the average FMNF is about 790.449 Hz.In summary, under the requirement that the optimal locations of bearings are obtained in one minute, after conducting the two-factor factorial experiments and analysis of variation, the optimal parameters of the GA are decided as ( z n p = (highlighted in Table 4), which can be applied to determine the optimal locations of bearings not only for the original design and the engineering changes, but also for the sensitivity studies of the design parameters.
In the 30 outcomes from the experiment of SDP-GA with optimal parameters ( ) ( ) , , 16, 74, 0.01 g m z n p = , the best and worst ones are 794.625Hz and 783.589Hz, which compare to the results obtained by the exhaustive search, the differences are about +0.000% and −388%, respectively.It is noticeable that the best result is slightly better than that of exhaustive search (794.622Hz) since the precision of the former is finer than the later, and therefore, we replace the value of 1 ω * with 794.625 Hz.Table 6 summarizes the optimal values of the four design variables, the corresponding FMNF's, and the final spindle-bearing system drawings for the exhaustive search and the best and worst runs.Moreover, since the average and standard deviation of FMNF for this combination are 790.449Hz and 3.749 Hz respectively, the approximate 95% confidence interval on the mean FMNF ( ) Therefore, if we can tolerate the potential errors of the SDP-GA approach as discussed above and indicated in Table 6, the number of times to run the subroutine for solving the eigenvalues can be reduced to 1200, which is equal to 6 2.893 10 − × of that the exhaustive search requires.Apparently, the quality of the optimal solutions shall be better if we are allowed to perform more calculations in the GA optimization.To investigate the influences of calculating times on the results of SDP-GA by using the optimal parameter values 16 z = and 0.01 Equation ( 19) can be used to decide the number g n needed to reach a certain average value of FMNF.For example, if we want the results with an average value 793 Hz, the required g n can be calculated and rounded to an integer, as 1536, which is equal to 20.5 minutes if each eigenvalue calculation takes 0.05 seconds, since ( )

Conclusion
This paper has developed a GA-based optimization approach to search for the optimal locations of bearings installed on a spindle shaft to maximize the FMNF.First, an FEM dynamic model of the spindle-bearing system is formulated, and by solving the eigenvalue problem derived from the equations of motion, the natural frequencies of the spindle system can be acquired.Next, the mathematical model is built, which includes the objective function to maximize FMNF and the constraints to limit the locations of the bearings with respect to the geometrical boundaries of the segments they located and the spacings between adjacent bearings.A customized GA optimization procedure, SDP-GA, is designed to accommodate the dependent characteristics of the constraints in the mathematical model.To verify the proposed SDP-GA optimization approach, a four-bearing installation optimization problem of a illustrative spindle system is investigated.The results show that the SDP-GA provides good convergence for the optimization searching process.By applying DOE and ANOVA, the optimal values of GA parameters are determined as 16 z = , , under the restriction that the number of executing the eigenvalue calculation subroutine is 1200.The best and worst outcomes for the optimal GA parameter are 794.625Hz and 783.589Hz, and for the best one, the locations for the four bearings are ( ) 189.04, 296.00, 558.00, 574.00 mm.Compared to the exhaustive search with a maximum FMNF as 794.622Hz, the 95% confidence interval of mean FMNF obtained by the SDP-GA with optimal parameter values is [ ] 789.107, 791.791Hz; however, the calculation time consumed by SDP-GA is about 6 2.893 10 − × of that of the exhaustive search.A linear regression equation is also derived to estimate necessary calculation efforts with respect to the specific quality of the optimization solution.From the results of this illustrative example, we can conclude that the proposed SDP-GA optimization approach is effective and efficient.
and for each element i : where AX ρ is the axial mass density, l is the length of the element, r is the radius of the element, E is Young's modulus, I is the moment of inertia, κ is the Timoshenko's shear correction factor, A is the cross-section area, and G is the shear modulus of elasticity.
the radial stiffness of bearing i , and j is the joint station which the bearing is acted on.The static radial stiffness V of an angular-contact ball bearing can be represented by the bearing stiffness equation as a function of axial preload ( )

1 −M K , and the roots of which are the eigenvalues of 1 − 1 p
M K .Expansion of Equation (8) by cofactors results in a ( ) 4 + th-degree polynomial equation in λ , i.e. the eigen- values of 1 − evaluated, we can obtain values of the global variables x by the global-local transition equation as indicated in Equation (

2 )
Extract the values of variables for each individual by reading b n sets of m binary digits and decoding into real numbers through the SDP.The fitness values, FMNF, for all individuals are evaluated by solving Equation (8), where the locations of bearings necessitate to be the joint stations of the finite element model, which is coded and tested antecedently.
we assign the global and

4 .
the origin point is located on the left end of the spindle shaft.The type I constraints for the other three design variables can be obtained in a similar way and the results are shown in Figure Type II constraints specify the least spacing between adjacent bearing pairs.As shown in Figure 5, the least distance between 1,1

Figure 7 .
Figure 7. Graph of average FMNF versus numbers of members for all values of m p .

Table 1 .
Dimensions and element numbers for the spindle shaft.

Table 2 .
The specifications of the bearings.

Table 3 .
The constraint equations and valid regions of the design variables (all units are in mm).

Table 4 .
Based on the information provided in

Table 4 ,
the ANOVA table for the two-factor problem is shown in

Table 4 .
Results of factorial design for parameters of SDP-GA.

Table 5 .
The analysis of variance table for the FMNF data.
as obtained in the above analysis, a series of 30-run tests are conducted, with total eigenvalue calculating times 3600 From the table, we can find that as the g n increases, the average FMNF increases while the standard deviation decreases, which confirms that the qualities of optimal solutions are improved with more iterations by the SDP-GA.The maximum and minimum FMNF obtained in each g n are also provided in the table, and they also perform better in larger g n .Finally, the last column of the table records the numbers of outcomes hitting 1 ω * , where we can find that the frequencies are greater with larger g

Table 7 .
Results of 30 runs for 74