Information feedback self-adaptive harmony search algorithm for the bovine cortical bone vibration-assisted drilling optimization

Measurement 149 (2020) 107020

Contents lists available at ScienceDirect

Measurement journal homepage: www.elsevier.com/locate/measurement

Information feedback self-adaptive harmony search algorithm for the bovine cortical bone vibration-assisted drilling optimization Li Shaomin, Zhang Deyuan ⇑, Shao Zhenyu, Tang Hui School of Mechanical Engineering and Automation, Beihang University, Beijing 100191, China Beijing Advanced Innovation Center for Biomedical Engineering, Beihang University, Beijing 100083, China

a r t i c l e

i n f o

Article history: Received 21 March 2019 Received in revised form 5 July 2019 Accepted 1 September 2019 Available online 5 September 2019 Keywords: Drilling optimization Harmony search Thrust force Bone drilling Vibration-assisted drilling

a b s t r a c t In orthopedics and surgery, bone drilling is an important operation. Increasing demands call for bone drilling force prediction and minimization. Vibration-assisted drilling (VAD) provides an opportunity to reduce the force. In the study, the information feedback self-adaptive harmony search (IFSHS) algorithm is presented for the thrust force minimization in bovine cortical bone VAD. The IFSHS improves the drilling parameters which leads to the optimal cutting state. The dynamic searching parameters, optimization parameters evaluation, and evaluation information feedback are applied to search the optimal solution. This force model was developed and employed with IFSHS. The VAD drilling experiments and parameters analysis are carried out. Experimental results have shown that the VAD has the potential to reduce the thrust force. Compared with conventional drilling (CD), the thrust force decreased by 18.4%–33.2% under suitable VAD parameters. Compared with the results of the original VAD parameters, the optimized thrust force decreased by 20.19%. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Bone cutting is performed around the world in the medical field which includes drilling, scraping, grooving, sawing [1]. The bone drilling operation is widely used in medical treatment which includes spine surgery, craniofacial surgery, fracture fixation, dental implantation and so on. In previous bone cutting studies, numerous attempts have been made to analyze the bone cutting process and reduce tissue damage. The cortical bone is a semibrittle material, and its anisotropic characteristics were considered for bone cutting force, chip formation mechanism, and temperature prediction [2,3]. A milling cutter with micro-cutting edges was designed to decrease surface damage [4]. Recently, prediction and optimization of the thrust force are concentrated on bone drilling studies, which is critical to a satisfactory result of the surgical [5]. Vibration-assisted drilling (VAD) is widely used in precision manufacturing, which could reduce cutting forces, temperature, exit burr height, and drilling-induced delamination of CFRP [6,7]. The large thrust force may lead to bone crack formation and damage to surrounding soft tissue, which could increase the recovery time for surgery. It is anticipated that VAD could drill holes with ⇑ Corresponding author at: School of Mechanical Engineering and Automation, Beihang University, Beijing 100191, China. E-mail address: [email protected] (D. Zhang). https://doi.org/10.1016/j.measurement.2019.107020 0263-2241/Ó 2019 Elsevier Ltd. All rights reserved.

the minimal thrust force to avoid unnecessary bone damage [8,9]. Alam investigated the VAD process of fresh bovine cortical bone [10]. This study showed that the force dropped significantly when vibration was applied. Wang compared the conventional drilling (CD) and low-frequency VAD in terms of bone microcracks [11]. The result showed that the micro-cracks produced during the VAD process were less and shorter than those produced by the CD. Alam assessed the impact of vibration-related parameters for the bone VAD [12]. The study indicated that vibration frequency has a greater impact on thrust force compared with vibration amplitude. The kinetic characteristic and parameters optimization of the VAD differ from the CD because of its complexity. The actual cutting angles of VAD are dynamic due to the superposition of the vibration and the feed motion. In the VAD process, the cutting parameters interact with each other. The relationship between the drilling parameters and drilling force is nonlinear. Under different constraints, there are different optimal parameter combinations in the machining process. Therefore, modeling and VAD parameters optimization of bone cutting is significant to accurate holes drilling without mechanical and thermal damage in the surgical operation. Analytical modeling, 3D finite-element modeling, and RSM method were applied to analyze the dynamic process prediction and optimization of bone drilling during the orthopedic surgery [13–15]. However, most of the previous bone drilling

2

S. Li et al. / Measurement 149 (2020) 107020

studies were the VAD and CD experimental analysis [16–18]. In these studies, few attempts were implemented to formulate the VAD thrust force model. Therefore, process analysis and modeling are important for bone VAD optimization. The optimization procedure has some challenges, such as multiple local optima, high computing complexity, and long computing time. Nowadays, meta-heuristic algorithms are proposed to face these challenges [19]. Harmony search (HS) is a novel intelligent search method [20]. It simulates the process of music creation for the optimization procedure [21]. In some of the previous works, the HS was applied to the multi-objective optimization problems for its excellent performance [22], such as the design optimization of mechanical arm control [23], max-cut problem [24], burr height minimization of VAD [7], multi-objective optimization of composite box beam [25], water network design [26], structural optimization [27], aircraft stiffened panels [19] and others [28,29]. Recent studies have shown that the Harmony Search algorithm has high efficiency and wide application. An improved HS was implanted in the Asymmetric Traveling Salesman Problem (ATSP), and its effectiveness of the technique was increased by almost 59% [30]. The harmony search algorithm still needs some improvement. The algorithm directive property exposed unsatisfactory during the local optimization, which decreases its efficiency [31,32]. The researchers proposed some modifications methods to improve search efficiency. These improvements were based on valuable information hidden in harmony memory [33], on-line variablefidelity surrogate-assisted mechanism [34], parallel chaotic local search enhanced mechanism [35], self-adaptive harmony particle swarm scheme [36], and efficient job sequence mapping with variable neighborhood search mechanism [37]. However, previous studies rarely consider the feedback of solution results, and the optimization process mainly relies on the search mechanism of the algorithm itself and the iterative updating of the harmony memory database. This paper proposes an information feedback self-adaptive harmony search algorithm (IFSHS) to improve the local optimization efficiency and optimize the bone VAD process. The present work is focused on the optimization of thrust force predictive model parameters using IFSHS algorithm in the bovine cortical bone VAD. In order to facilitate computations, the proposed thrust force model is simplified into the primary cutting edges force, and its accuracy improved by adopting an empirical coefficient. Some standard benchmark functions are used to find out the favorable algorithm parameters ranges. For the thrust force minimization, the VAD force model is combined with IFSHS method to search the optimal drilling parameters. Finally, the optimization experiments are carried out, and the IFSHS is applied to the cutting process improvement to minimize the thrust force. 2. The IFSHS algorithm HS algorithm is an intelligent optimization method. Its principle is derived from music creative practice. In this optimization procedure, every N-dimensional optimal solution is regarded as one beautiful harmony, and this optimization procedure is regarded

as the musician’s improvisation. Fig. 1 shows a harmony creative process from some particular value ranges. The HS process for the optimal solution is as follows [20,21]: Step 1. Define the engineering problem and the HS parameters. Step 2. Harmony Memory (HM) assignment. Step 3. Generate a new harmony. Step 4. Generate new HM. Step 5. Inspection termination condition. The HMCR and PAR influence the global and local searching efficiency, respectively. The PAR and BW affect the convergence speed in the multi-objective optimization problem. Therefore, the reasonable parameter configuration is crucial to the efficient operation of the algorithm. In conventional HS, its parameters are constant values. To increase the searching efficiency of the HS, Mahdavi et al. presented a novel improvement idea which is variable algorithm parameters [31]. This method called the improved HS algorithm (IHS). Partial parameters of IHS vary with the search number for the improvement of the convergence efficiency. They could be given as follows:

8
 < PARðtÞ ¼ PARmin þ ðPARmax  PARmin Þ NIt 
   BWmin t : BWðtÞ ¼ BWmax exp Ln BW NI max

where PAR(t) is the PAR value for the iterative variable t, PARmax and PARmin are the available extreme value of PAR, BW(t) is the bandwidth value for the iterative variable t, BWmax and BWmin are the available extreme value of bandwidth, and NI is the maximum search number. In the IHS, PAR and BW vary according to the given rules. However, the solution evaluation information according to the objective function was not considered. In addition, the fine tuning of HMCR is not considered in the IHS, which could avoid the local optimum. In order to improve these problems, the parameters are given as follows:

h 8
t  HMCR i > HMCRðtÞ ¼ HMCR exp nðtÞ þ Ln HMCRmax > min NI > min <
t  PARðtÞ ¼ PARmin þ ðPARmax  PARmin Þ NI þ nðtÞ > h   i > > : BWðtÞ ¼ BW exp nðtÞ þ
t Ln BWmin max BWmax NI

ð2Þ

where HMCR(t) is the HMCR value for the iterative variable t, HMCRmax and HMCRmin are the available extreme value of the HMCR, nðtÞ is the feedback factor for the iterative variable t, the feedback factor nðtÞ records the evaluation information for the iterative variable t, and the nðtÞ is given as follows:

8 ¼0 < t < 3 ) nðtÞ  nðtÞ ¼ expðjgðt  1ÞjÞgðt  1Þ > 0 :t > 3 ) nðtÞ ¼ expðjgðt  1ÞjÞgðt  1Þ  0

gðt  1Þ ¼

Fig. 1. The HS algorithm principle [7].

ð1Þ

rðt  1Þ  rðt  2Þ rðt  2Þ

ð3Þ

ð4Þ

where g(t-1) is the search result variation ratio for the iterative variable t-1, and r(t-1) is the search result for the iterative variable t-1. The IFSHS algorithm employs a self-adaptive parameters adjustment program which includes dynamic parameters and feedback factors. Fig. 2 shows the optimal solution search process of the IFSHS algorithm.

3

S. Li et al. / Measurement 149 (2020) 107020

Fig. 2. The search process of the IFSHS algorithm.

3. Bone VAD force model Bone cutting is an important operation in medical treatment. In bone drilling, the drill bit includes two primary cutting edges and one chisel edge [38,39]. They are the source of cutting force. According to previous experimental data [40], the total thrust force is approximately twice as large as the force of the primary cutting edge (Fig. 3). Therefore, a simplified force model is given as follows:

F th ¼ F th1 þ F th2  2F th1

ð5Þ

where Fth, Fth1, and Fth2 are the total force, primary cutting edge force, and the chisel edge force, respectively. The orthogonal cutting model could be used for the calculation of the thrust force [41,42]. The thrust force of bone drilling could be calculated as follows [7]:

F th1 ¼ 2

N X

DF th1

ð6Þ

DF t1 ¼

ks hd DLsinðknd  cnd Þ sin/nd cosð/nd þ knd  cnd Þ

ð8Þ

DF c1 ¼

ks hd DLcosðknd  cnd Þ sin/nd cosð/nd þ knd  cnd Þ

ð9Þ

gd ¼ tan1 DL ¼

ð7Þ

ð10Þ

Dcosðsin1 xÞ  2Tcotðp  wÞ 2Nsinp

ð11Þ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 1 DL sin p þ T 2 ri ¼ ½Rcosðsin1 xÞ  i  2

x¼

ð12Þ

T R

i¼1

DF th1 ¼ DF t1 sin pcos gd þ DF c1 sin gd

f 2pr i

cn ¼ tan1

ð13Þ " 2 ðq2i  x2 sin pÞtan b0

Fig. 3. Cutting edge force model.

ðq2i  x2 Þ

1=2

sin p



xcos p 1=2 ðq2i  x2 Þ

# ð14Þ

4

qi ¼

S. Li et al. / Measurement 149 (2020) 107020

ri R

i ¼ 1; 2; 3  

ð15Þ

p cnd

knd ¼

þ

6

/nd ¼

ð16Þ

2

p cnd 6

þ

ð17Þ

4

cnd ¼ cn þ gd

ð18Þ

hd ¼ hf sinp=cosgd

ð19Þ

where DFth1 is the elemental thrust force; ks is the shear strength; hd is the dynamic uncut chip thickness; DL is the element width; ri is the element radius; f is the feed rate; /nd is the dynamic shear angle; gd is the dynamic feed angle; kn and knd are the normal and dynamic friction angle, respectively; cn and cnd are the normal and dynamic rake angle, respectively; 2p is the point angle; w is the chisel edge angle; b0 is the helix angle; R is the drill radius; and 2 T is the web thickness (Fig. 4). The VAD model is showed in Fig. 1. Dynamic cutting geometry calculation is an essential problem for the thrust force evaluation. The dynamic drill bit axial displacement and the dynamic uncut chip thickness hf due to the vibration are focused. The instantaneous position x(t) and is given by [43,44]:

xðtÞ ¼ Asinð2pktÞ þ fnt

ð20Þ

where A and k are the vibration amplitude and frequency, respectively; and n is the spindle speed (Fig. 5).

h ¼ 2pnt

ð21Þ 

xðhÞ ¼ Asin

k f h þ h n 2p

ð22Þ ð23Þ

where m is the minimum which satisfies the following condition:

In order to confirm the value ranges of IFSHS parameters and test algorithm performance, some experiments are implemented based on various standard benchmark problems. They are Sphere function, Rosenbrock function, Step function, Rastrigin function, Ackley function, Rotated hyper-ellipsoid function, and Six-Hump Camel-Back function. The details are given as follows: A. Sphere function (Fig. 6(a)), and it could be expressed as:

xðhÞ  xmax ðh  mpÞ if

min f ðxÞ ¼

n X

x2i

 100  x  100

ð26Þ

The optimal solution is x = 0, f (x) = 0. B. Rosenbrock function (Fig. 6(b)), and it could be expressed as:

ð24Þ

In VAD, the uncut chip thickness is a dynamic value. It is related to material removal efficiency and thrust force calculation. If the axial position of the drill bit is under the machined surface, it is an effective cutting. Otherwise, the tool does not cut the bone. It could be expressed as follows:

0

4.1. Application to standard test problems

i¼1

xðh - mpÞ > xðh  ðm þ 1ÞpÞ

hf ¼

4. Bovine cortical bone VAD optimization



xmax ðhÞ ¼ xðh  mpÞ



Fig. 5. Kinematics of vibration-assisted drilling.

xðhÞ > xmax ðh  mpÞ otherwise

ð25Þ

min f ðxÞ ¼

n1 X 2 ½100ðxiþ1  x2i Þ þ ðxi  1Þ2 

 30  x  30

ð27Þ

i¼1

The optimal solution is x = (1, 1,. . ., 1), f (x) = 0. C. Step function (Fig. 6(c)), and it could be expressed as:

min f ðxÞ ¼

n X i¼1

Fig. 4. The geometry of a twist drill bit.

ðbxi þ 0:5cÞ2

 100  x  100

ð28Þ

5

S. Li et al. / Measurement 149 (2020) 107020

Fig. 6. The standard benchmark problems in the experiments.

The optimal solution is x = 0, f (x) = 0. D. Rastrigin function (Fig. 6(d)), and it could be expressed as:

min f ðxÞ ¼

n X ½x2i  10cosð2pxi Þ þ 10 i¼1

The optimal solution is x = 0, f (x) = 0.

 5:12  x  5:12

ð29Þ

6

S. Li et al. / Measurement 149 (2020) 107020

Table 1 The result under different HMS values (n = 5).

A B C D E F G

Optimal solution

HMS 5

10

30

50

0 0 0 0 0 0 1.0316285

0.000000 ± 0.000000 0.178227 ± 0.038329 0.000000 ± 0.000000 0.000000 ± 0.000000 0.000000 ± 0.000000 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000087 ± 0.000138 0.335447 ± 0.185159 0.100000 ± 0.307793 0.018564 ± 0.009088 0.019509 ± 0.005182 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000139 ± 0.000056 0.459206 ± 0.108291 0.200000 ± 0.410391 0.018700 ± 0.008832 0.019655 ± 0.005234 0.000006 ± 0.000016 1.031628 ± 0.000000

0.000167 ± 0.000069 0.757607 ± 1.112136 0.250000 ± 0.444261 0.021484 ± 0.008987 0.021200 ± 0.005164 0.000005 ± 0.000010 1.031628 ± 0.000000

0.1–0.9

0.3–0.9

0.5–0.9

0.7–0.9

0.000160 ± 0.000055 3.863342 ± 4.613915 0.100000 ± 0.307793 0.022383 ± 0.007901 0.022898 ± 0.004456 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000149 ± 0.000060 11.41654 ± 37.38148 0.200000 ± 0.410391 0.025796 ± 0.010180 0.022065 ± 0.004745 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000139 ± 0.000050 2.271712 ± 4.249604 0.150000 ± 0.366347 0.027447 ± 0.007556 0.022459 ± 0.005737 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000128 ± 0.000069 4.669142 ± 5.934000 0.130000 ± 0.550119 0.022759 ± 0.009694 0.021648 ± 0.006193 0.000000 ± 0.000000 1.031628 ± 0.000000

0.1–0.9

0.3–0.9

0.5–0.9

0.7–0.9

0.000155 ± 0.000047 4.053078 ± 6.064303 0.200000 ± 0.410391 0.023700 ± 0.013086 0.023222 ± 0.005669 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000131 ± 0.000059 13.11580 ± 36.10206 0.250000 ± 0.444261 0.026931 ± 0.014512 0.020495 ± 0.004611 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000129 ± 0.000049 4.001159 ± 4.630073 0.350000 ± 0.670820 0.027812 ± 0.011666 0.020711 ± 0.005866 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000166 ± 0.000042 3.858394 ± 5.211834 0.300000 ± 0.470162 0.024636 ± 0.009727 0.021397 ± 0.005857 0.000000 ± 0.000000 1.031628 ± 0.000000

Table 2 The result under different HMCR values (n = 5). Optimal solution

A B C D E F G

0 0 0 0 0 0 1.0316285

HMCR

Table 3 The result under different PAR values (n = 5). Optimal solution

A B C D E F G

0 0 0 0 0 0 1.0316285

PAR

Table 4 The result under different BW values (n = 5).

A B C D E F G

Optimal solution

BW 0.001–0.5

0.005–0.5

0.01–0.5

0.05–0.5

0.1–0.5

0 0 0 0 0 0 1.0316285

0.000000 ± 0.000000 5.982693 ± 6.514515 0.150000 ± 0.366347 0.000009 ± 0.000004 0.000437 ± 0.000119 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000001 ± 0.000000 3.853278 ± 5.105011 0.200000 ± 0.410391 0.000246 ± 0.000091 0.002286 ± 0.000501 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000007 ± 0.000003 3.703102 ± 6.142217 0.100000 ± 0.307793 0.001014 ± 0.000514 0.004294 ± 0.001310 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000154 ± 0.000058 4.669142 ± 5.934000 0.050000 ± 0.223606 0.027424 ± 0.009601 0.022586 ± 0.005533 0.000000 ± 0.000000 1.031627 ± 0.000001

0.000462 ± 0.000217 12.093894 ± 38.728858 0.200000 ± 0.523148 0.148804 ± 0.220681 0.048033 ± 0.011299 0.000000 ± 0.000000 1.031626 ± 0.000002

Table 5 The result under different BW values (n = 5).

A B C D E F G

Optimal solution

BW 0.001–0.1

0.001–0.3

0.001–0.5

0.001–0.7

0.001–0.9

0.001–1.0

0 0 0 0 0 0 1.0316285

0.000000 ± 0.000000 3.573126 ± 5.161816 0.100000 ± 0.307793 0.000011 ± 0.000005 0.000487 ± 0.000110 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000000 ± 0.000000 9.371705 ± 24.00213 0.150000 ± 0.366347 0.000011 ± 0.000005 0.000437 ± 0.000115 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000000 ± 0.000000 2.516559 ± 2.742779 0.100000 ± 0.307793 0.000010 ± 0.000004 0.000504 ± 0.000103 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000000 ± 0.000000 2.279634 ± 2.823624 0.300000 ± 0.470162 0.000010 ± 0.000004 0.000449 ± 0.000116 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000000 ± 0.000000 4.247538 ± 4.735386 0.200000 ± 0.410391 0.000011 ± 0.000003 0.000474 ± 0.000107 0.000000 ± 0.000000 1.031628 ± 0.000000

0.000000 ± 0.000000 2.913955 ± 3.399329 0.150000 ± 0.366347 0.000012 ± 0.000003 0.000437 ± 0.000119 0.000000 ± 0.000000 1.031628 ± 0.000000

7

S. Li et al. / Measurement 149 (2020) 107020

E. Ackley function (Fig. 6(e)), and it could be expressed as:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 ! u n n u1 X 1X min f ðxÞ ¼ 20exp@0:2t x2i A  exp cosð2pxi Þ n i¼1 n i¼1 0

þ 20 þ e

 32  x  32

The optimal solution is x = 0, f (x) = 0.

ð30Þ

F. Rotated hyper-ellipsoid function (Fig. 6(f)), and it could be expressed as:

min f ðxÞ ¼

n i X X i¼1

!2 xj

 100  x  100

j¼1

The optimal solution is x = 0, f (x) = 0.

Fig. 7. The standardized search results of the algorithm.

ð31Þ

8

S. Li et al. / Measurement 149 (2020) 107020

G. Six-Hump Camel-Back function (Fig. 6(g)), and it could be expressed as:

1 min f ðxÞ ¼ 4x21  2:1x41 þ x61 þ x1 x2  4x22 þ 4x42 3

5x5 ð32Þ

The optimal solution is x = (0.08983, 0.7126), f (x) = 1.0316285. Some of these functions are the single-peaked function, and some are multimodal functions. The search results are the average value of 20 simulations. The search times is 50,000. The logical value ranges of the algorithm parameters are investigated, as shown in Tables 1–5. In order to analyze the calculation results under different parameters, and the search results normalization were conducted (Fig. 7). The normalization formula is as follows:

Ni ¼

ai  Amin Amax  Amin

methods. The search results of Rosenbrock is poor, but the result is still better than other algorithms. According to the above analysis, the applicable value ranges of the algorithm parameters are summarized in Table 6. Based on the above parameter setting, the convergence curves for the standard benchmark optimization problems are shown in Fig. 8. The figure demonstrates that the IFSHS algorithm could reach the GO after 20,000 iterations. 4.2. Force optimization of bovine cortical bone VAD Bovine cortical bone VAD experiments were implemented on a VAD center. It includes a drilling machine, a driver controller, computer, and axial vibrator (Fig. 9). In the experiments, the vibration

ð33Þ

where Ni is the normalized value, ai is the initial value, Amax and Amin are the maximum and minimum values in the initial results, respectively. Fig. 7(a) shows that the higher HMS leads to an increase in the optimal value, which indicates that large HMS is not conducive to the search of the optimal solution. Change in the HMCR value range affect the search results, but the influence does not show linear characteristics. The larger value (0.5–0.9) may lead to a decrease in the search results. The effect of HMCR and PAR on the results show similar characteristics. The left boundary of BW value has a more significant influence on the search results (Fig. 9(d) and (e)). The small left boundary value is conducive to searching for the optimal solution. The optimum value of parameters exists distinction for different problems. However, the logical value ranges could be acquired. Small size HMS (5–10) is recommended, which could decrease the calculated amount. As shown in Tables 2 and 3, larger HMCR and PAR could improve the optimization performance. According to the experiment results, the IFSHS gains a stable performance when the HMCR and PAR are between 0.5 and 0.9. At the initial stage, large BW could increase the diversity of HM. In the final stage of the search, small BW could contribute to convergence speed [31], as shown in Tables 4 and 5. The BW value range chooses 0.001 to 0.5. The comparison results of HS, GA, and IFSHS are shown in Table 7, which is used to analyze the IFSHS algorithm performance. The dimension of the functions is 30. Table 7 shows that IFSHS improves the optimum solution significant compared to other

Fig. 8. Typical convergence curves for the standard benchmark optimization problems.

Table 6 Parameters of IFSHS. Parameter

HMS

HMCR

PAR

BW

Value range

10

0.5–0.99

0.5–0.99

0.001–0.5

Fig. 9. Bovine femur bone VAD experiment platform.

Table 7 Optimization results of HS, GA, and IFSHS for benchmark functions [45,46].

A B C D E F

Optimal solution

HS

GA

IFSHS

0 0 0 0 0 0

2.81E04 ± 1.25E05 3.29E+01 ± 2.88E+01 0.00E+00 ± 0.00E+00 4.66E02 ± 4.92E03 1.24E02 ± 9.94E05 4.30E+03 ± 1.36E+03

2.13E04 ± 1.07E05 1.35E+01 ± 1.24E+01 0.00E+00 ± 0.00E+00 5.32E02 ± 3.54E03 3.46E02 ± 2.65E05 3.65E+03 ± 2.52E+03

2.13E04 ± 1.07E05 1.25E+01 ± 1.52E+01 0.01E02 ± 0.02E02 2.54E02 ± 1.43E03 1.01E02 ± 1.98E05 2.74E+03 ± 1.65E+03

S. Li et al. / Measurement 149 (2020) 107020

9

Fig. 10. Thrust force of various drilling conditions. (a): Spindle speed 4000 r/min, vibration frequency 200 Hz; (b): Vibration amplitude 0.002 mm, spindle speed 4000 r/min; (c): Spindle speed 4000 r/min, feed rate 0.025 mm/r; (d): Feed rate 0.025 mm/r, vibration amplitude 0.002 mm.

Table 8 The drilling parameters value range. Parameters Values

f mm/r

n r/min

k Hz

A

0.01–0.02

1500–3000

100–1000

2.3–3.5

lm

amplitude is optional between 2.3 lm and 3.5 lm, and the frequency is adjustable from 0 to 1000 Hz. Fresh bovine femurs were used in this bone drilling experiment, which was collected within 24 h after slaughter. All experiments used the same batch of bovine femurs. The central parts of the femur were used as the test sample, and the sample wideness was 50 mm. The thrust force of the drilling experiment was measured and collected in the direction perpendicular to the axis of the bovine femur. The force was

Fig. 11. Experimental values for thrust force.

10

S. Li et al. / Measurement 149 (2020) 107020

Table 9 The bone drilling experiments. Parameters

f mm/r

n r/min

k Hz

A

Values

0.01/0.02

1500/2000/3000

0/200/400/600/800/1000

2.3/3.5

Fig. 12. Experimental values for thrust force.

lm

S. Li et al. / Measurement 149 (2020) 107020

measured and recorded by a dynamometer system. It included a three-dimensional dynamometer (Kistler 9254), a multichannel charge amplifier (Type 5070), and a data acquisition system (Type 5697). The resolution of the system was 0.02 N, and its sampling rate was 1000 Hz. The 3D surface plots were applied to investigate the VAD process parameter selection during bone drilling. Fig. 10 indicates the thrust force response under the various combination of the drilling parameters. Fig. 10(a) shows a positive and nearly linear correlation between feed rate and thrust force. A similar tendency could be seen in Fig. 10(b). The result is reasonable because a larger feed rate lead to a higher material removal rate, and it could increase the cutting force. In addition, this process could also result in an increase in cutting temperature, which deteriorates the cutting conditions. Vibration machining significantly reduces the thrust force, but the force is not sensitive to changes in vibration amplitude. With the increase of amplitude, the thrust force increases slightly. Fig. 10(c) shows that the interaction between vibration frequency and amplitude is obvious. At higher vibration amplitude, the thrust force is reduced when the vibration frequency is within the intermediate value range. Fig. 10(d) indicates that higher spindle speed could help reduce the thrust force. High spindle speed may lead to an increase of contact surfaces friction. However, it could decrease the feed rate which results in a lower thrust force. According to the above results, low feed rate, low amplitude, and high drilling speed are recommended. Optimal vibration frequency under different conditions may require targeted analysis. Although the influence of various parameters on the force has a certain regularity, it still shows non-linear characteristics when multiple parameters work together. In VAD, the optimum vibration condition (vibration frequency and amplitude) under different cutting condition are various. Besides, minimum force may not appear at the achievable boundary value of vibration frequency. These may be due to the nonlinear characteristic. Therefore, the multiobjective optimization tool is required for the thrust force minimization. The IFSHS and force model was applied to improve the bovine cortical bone VAD process. The optimization procedure tries to find suitable parameters under certain tool geometry characteristics, which could minimize the thrust force. The algorithm running process is shown in Fig. 2. Its termination condition is that the number of cycles reaches 20,000. Table 8 shows restrictions of drilling parameters. The HM is initialized with uniform interval values within the value range. The HSS twist drill was applied to the bovine femur cortical bone drilling experiments. Its diameter was 2.4 mm with 125° point angle and 25° helix angle. The chisel edge angle was 135° and the chisel edge length was 0.5 mm. The optimal drilling parameters from the IFSHS optimization experiment are spindle speed 3000 r/min, feed rate 0.01 mm/r, vibration frequency 236 Hz, and vibration amplitude 3.5 lm. The thrust force is 10.67 N. According to the results of the optimization parameters, Fig. 11 shows the experimental thrust force and its reduction percentage at different vibration frequencies. Fig. 11(a) and (b) show the results for 2.3 lm and 3.5 lm amplitude, respectively. The results show that the optimization experiment is consistent with the cutting test. Low feed rate and high drilling speed could contribute to the reduction of the thrust force. The recommended vibration amplitude and frequency of simulation are 3.5 lm and 236 Hz, respectively. Suitable vibration parameters improve the separation stability of the cutting edge and the materials, which could reduce the drilling heat during the bone cutting. Dynamic uncut chip thickness and intermittent cutting could reduce the force. The deviation of the optimal result and the experiment value is within

11

7.63%. The VAD reduces the thrust force by 18.4%–33.2% compared to CD. The thrust force of 200 Hz and 1000 Hz are approximate. This shows that a higher vibration frequency may result in a higher thrust force reduction. It is consistent with the model analysis that the higher vibration frequency could result in larger thrust force decrease. Comparing the results of the original VAD and the optimized VAD, the thrust force value decreased 20.19%. Some bone drilling experiments were carried out to analyze the search results. The experimental parameters are shown in Table 9. The statistical analysis of the experiments and the prediction of the optimal solution are shown in Fig. 12. Fig. 12(a) and (b) indicates that the minimum thrust force may occur under the optimal vibration frequency condition (200–300 Hz). This indicates that the optimal parameter may not appear at the boundary position of the value range. The higher spindle speed could result in a decrease in the thrust force (Fig. 12(b)). The feed rate is reduced at higher spindle speed (3000 r/min), which could lead to a lower material removal rate. The statistical results also show that a higher feed rate (0.02 mm/r) could result in a higher thrust force, while large vibration amplitude (3.5 lm) could result in a decrease of force. It’s because the higher feed rate could increase the material removal rate. However, the large vibration amplitude may improve the drilling condition due to the higher separation rate of tool and material. The results of statistical analysis are consistent with the search results of the algorithm. 5. Conclusion The application of IFSHS algorithm optimization for the thrust force reduction in VAD is presented in this paper. The VAD process parameters were considered for the thrust force mathematical models development. IFSHS implemented parameters optimization for specific drill bit diameters and obtained the minimum thrust force. Some important modifications are applied to the parameters optimization in the bone drilling process. In order to simplify the calculation, the proposed thrust force model is simplified into the primary cutting edges force, and its accuracy is improved by adopting an empirical coefficient. The feedback factor nðtÞ is employed to the fluctuation mechanism of algorithm parameters, and the fluctuation mechanism is applicable to HMCR, PAR, and BW. VAD parameters interact with each other, and the effects of these parameters on the target value are nonlinear. These characteristics indicate the necessity of using intelligent algorithms for the VAD parameters optimization. The dynamic searching parameters, optimization parameters evaluation, evaluation information feedback, and diverse harmony memory are applied to search the optimal VAD machining parameters in the bone drilling. The results indicated that the models are effective in predicting and minimizing thrust force in VAD. The results show that the vibration amplitude and feed rate are important factors for the VAD thrust force minimization. The optimal vibration conditions are different for various drilling parameters combinations. This study also reveals that a combination of appropriate vibration frequency and amplitude, low feed, and high spindle speed are essential in controlling thrust force. Appendix The detailed process of the HS algorithm is described below [20]: Step 1. Define the objective function and HS parameters. This engineering problem is analyzed, and the objective function is obtained. This objective function which could be formulated as:

12

S. Li et al. / Measurement 149 (2020) 107020

minimizing f ðxÞ subject to

xi 2 X i ;

i ¼ 1; 2; . . . ; N

ð34Þ

where f (x) is a particular object function, x is the independent variable, N is the solution vector dimension, and Xi is the achievable range of the independent variable xi. In order to implement this search operation, this step assigns values to some algorithm parameters based on the engineering problem characteristic. The Harmony Memory Size (HMS) is the number of the optimal solution set. Number of Improvisations (NI) is the number of searches. Harmony Memory Considering Rate (HMCR) and Pitch Adjusting Rate (PAR) are two probability parameters for searching optimal solution. The Band Width (BW) is the step size parameter to adjust the optimal solution. The parameters configuration is crucial to search performance. Therefore, the logical parameter value range should be analyzed. Step 2. Harmony Memory (HM) assignment. The HM includes various harmony vector, and its capacity is HMS. Every harmony vector (x1, x2,. . ., xN) represents one solution to the problem. The diversity of HM is significant for the searching efficiency of the optimal solution. The HM expression could be given as follows:

2

x11 6 2 6 x1 6 HM ¼ 6 . 6 .. 4 xHMS 1

x12 x22 .. . xHMS 2

... ... ... ...

x1N x2N .. .

3 7 7 7 7 7 5

ð35Þ

xHMS N

Step 3. Generate a new harmony. The improvisation is used to create a new vector (x1, x2,. . ., xN). This process adopts three methods. One is the memory considerations, its rule is that the HS chooses a harmony from the current HM. The second is the pitch adjustments which means that the HS slightly adjusts the vector value from the HM. The third is the random selection which denotes that the HS randomly chooses the harmony within a certain range. The procedure could be summarized as follows:

( xnew i

 2 x1i ; x2i ; . . . ; xHMS xnew i i 2 Xi xnew i

with probability HMCR; with probability 1-HMCR: ð36Þ xnew i

Every independent parameter chosen from the HM may be adjusted, and the probability is PAR. The procedure could be expressed as follows:

 xnew i

¼ xnew  r 2  BW xnew i i new xi

with probability PAR; with probability 1-PAR:

ð37Þ

where r2 is the random number, and its value range is 0 to 1. Step 4. Generate new HM. In the search procedure, the algorithm estimates the current solutions in the HM and the new one. These better harmonies are left in the HM, and the worse one is deleted. These cycles constantly improve the harmony quality of HM. Step 5. Inspection termination condition. If the searching result satisfies the termination condition, then stop running. Otherwise, the algorithm repeats step 3 and step 4. References [1] R.K. Pandey, S.S. Panda, Multi-performance optimization of bone drilling using Taguchi method based on membership function, Measurement 59 (2015) 9– 13. [2] Z. Liao, D. Axinte, D. Gao, On modelling of cutting force and temperature in bone milling, J. Mater. Process. Technol. 266 (2019) 627–638. [3] Z. Liao, D.A. Axinte, On chip formation mechanism in orthogonal cutting of bone, Int. J. Mach. Tools Manuf 102 (2016) 41–55.

[4] Z. Liao, D.A. Axinte, D. Gao, A novel cutting tool design to avoid surface damage in bone machining, Int. J. Mach. Tools Manuf 116 (2017) 52–59. [5] R.K. Pandey, S.S. Panda, Optimization of bone drilling parameters using greybased fuzzy algorithm, Measurement 47 (2014) 386–392. [6] D. Geng, Y. Liu, Z. Shao, Z. Lu, J. Cai, X. Li, X. Jiang, D. Zhang, Delamination formation, evaluation and suppression during drilling of composite laminates: a review, Compos. Struct. 216 (2019) 168–186. [7] L. Shaomin, Z. Deyuan, G. Daxi, S. Zhenyu, T. Hui, Modeling and drilling parameters optimization on burr height using harmony search algorithm in low-frequency vibration-assisted drilling, Int. J. Adv. Manuf. Technol. 101 (9– 12) (2018) 2313–2325. [8] Z. Li, D. Yang, W. Hao, T. Wu, S. Wu, X. Li, A novel technique for micro-hole forming on skull with the assistance of ultrasonic vibration, J. Mech. Behav. Biomed. Mater. 57 (2016) 1–13. [9] Z. Li, D. Yang, W. Hao, S. Wu, Y. Ye, Z. Chen, X. Li, Ultrasonic vibration-assisted micro-hole forming on skull, Proc. Instit. Mech. Eng., Part B: J. Eng. Manuf. 231 (14) (2015) 2447–2457. [10] K. Alam, A.V. Mitrofanov, V.V. Silberschmidt, Experimental investigations of forces and torque in conventional and ultrasonically-assisted drilling of cortical bone, Med. Eng. Phys. 33 (2) (2011) 234–239. [11] Y. Wang, M. Cao, Y. Zhao, G. Zhou, W. Liu, D. Li, Experimental investigations on microcracks in vibrational and conventional drilling of cortical bone, J. Nanomater. 2013 (2013) 1–5. [12] K. Alam, M. Ghodsi, A. Al-Shabibi, V. Silberschmidt, Experimental study on the effect of point angle on force and temperature in ultrasonically assisted bone drilling, J. Med. Biol. Eng. 38 (2) (2017) 236–243. [13] L. Qi, X. Wang, M.Q. Meng, 3D finite element modeling and analysis of dynamic force in bone drilling for orthopedic surgery, Int. J. Numer. Method Biomed. Eng. 30 (9) (2014) 845–856. [14] J. Sui, N. Sugita, K. Ishii, K. Harada, M. Mitsuishi, Mechanistic modeling of bonedrilling process with experimental validation, J. Mater. Process. Technol. 214 (4) (2014) 1018–1026. [15] S.R.H. Davidson, Drilling in bone: modeling heat generation and temperature distribution, J. Biomech. Eng. 125 (3) (2003) 305. [16] K.I.A.-L. Al-Abdullah, H. Abdi, C.P. Lim, W.A. Yassin, Force and temperature modelling of bone milling using artificial neural networks, Measurement 116 (2018) 25–37. [17] T. MacAvelia, A. Ghasempoor, F. Janabi-Sharifi, Force and torque modelling of drilling simulation for orthopaedic surgery, Comput. Methods Biomech. Biomed. Eng. 17 (12) (2014) 1285–1294. [18] T. Rajmohan, K. Palanikumar, S. Prakash, Grey-fuzzy algorithm to optimise machining parameters in drilling of hybrid metal matrix composites, Compos. B Eng. 50 (2013) 297–308. [19] B. Keshtegar, P. Hao, Y. Wang, Q. Hu, An adaptive response surface method and Gaussian global-best harmony search algorithm for optimization of aircraft stiffened panels, Appl. Soft Comput. 66 (2018) 196–207. [20] W.G. Zong, J.H. Kim, G.V. Loganathan, A new heuristic optimization algorithm: harmony search, Simul. Trans. Soc. Model. Simul. Int. 76 (2) (2001) 60–68. [21] K.S. Lee, Z.W. Geem, A new meta-heuristic algorithm for continuous engineering optimization: harmony search theory and practice, Comput. Methods Appl. Mech. Eng. 194 (36–38) (2005) 3902–3933. [22] M. Shabani, S. Abolghasem Mirroshandel, H. Asheri, Selective refining harmony search: a new optimization algorithm, Expert Syst. Appl. 81 (2017) 423–443. [23] Z. He, B. Pan, Z. Liu, X. Tang, The mechanical arm control based on harmony search genetic algorithm, Cluster Computing 20 (4) (2017) 3251–3261. [24] Y.-H. Kim, Y. Yoon, Z.W. Geem, A comparison study of harmony search and genetic algorithm for the max-cut problem, Swarm Evol. Comput. 44 (2019) 130–135. [25] K. Lakshmi, A.R.M. Rao, Multi-objective optimal design of composite box beam using hybrid adaptive harmony search with dynamically reconfigurable harmony memory, Proc. Instit. Mech. Eng., Part C: J. Mech. Eng. Sci. 233 (2) (2019) 713–734. [26] Z.W. Geem, Particle-swarm harmony search for water network design, Eng. Optim. 41 (4) (2009) 297–311. [27] M.-Y. Cheng, D. Prayogo, Y.-W. Wu, M.M. Lukito, A Hybrid Harmony Search algorithm for discrete sizing optimization of truss structure, Autom. Constr. 69 (2016) 21–33. [28] K. Abhishek, S. Datta, S.S. Mahapatra, Multi-objective optimization in drilling of CFRP (polyester) composites: Application of a fuzzy embedded harmony search (HS) algorithm, Measurement 77 (2016) 222–239. [29] S. Chatterjee, S.S. Mahapatra, K. Abhishek, Simulation and optimization of machining parameters in drilling of titanium alloys, Simul. Model. Pract. Theory 62 (2016) 31–48. [30] U. Boryczka, K. Szwarc, The harmony search algorithm with additional improvement of harmony memory for asymmetric traveling salesman problem, Expert Syst. Appl. 122 (2019) 43–53. [31] M. Mahdavi, M. Fesanghary, E. Damangir, An improved harmony search algorithm for solving optimization problems, Appl. Math. Comput. 188 (2) (2007) 1567–1579. [32] Q.-K. Pan, P.N. Suganthan, M.F. Tasgetiren, J.J. Liang, A self-adaptive global best harmony search algorithm for continuous optimization problems, Appl. Math. Comput. 216 (3) (2010) 830–848. [33] K. Luo, J. Ma, Q. Zhao, Enhanced self-adaptive global-best harmony search without any extra statistic and external archive, Inf. Sci. 482 (2019) 228–247.

S. Li et al. / Measurement 149 (2020) 107020 [34] J. Yi, L. Gao, X. Li, C.A. Shoemaker, C. Lu, An on-line variable-fidelity surrogateassisted harmony search algorithm with multi-level screening strategy for expensive engineering design optimization, Knowl.-Based Syst. 170 (2019) 1– 19. [35] J. Yi, X. Li, C.-H. Chu, L. Gao, Parallel chaotic local search enhanced harmony search algorithm for engineering design optimization, J. Intell. Manuf. 30 (1) (2019) 405–428. [36] F. Zhao, Y. Liu, C. Zhang, J. Wang, A self-adaptive harmony PSO search algorithm and its performance analysis, Expert Syst. Appl. 42 (21) (2015) 7436–7455. [37] F. Zhao, Y. Liu, Y. Zhang, W. Ma, C. Zhang, A hybrid harmony search algorithm with efficient job sequence scheme and variable neighborhood search for the permutation flow shop scheduling problems, Eng. Appl. Artif. Intell. 65 (2017) 178–199. [38] V. Chandrasekharan, S.G. Kapoor, R.E. Devor, A mechanistic approach to predicting the cutting forces in drilling: with application to fiber-reinforced composite materials, J. Eng. Ind. 117 (4) (1995) 559–570.

13

[39] A. Paul, S.G. Kapoor, R.E. DeVor, A chisel edge model for arbitrary drill point geometry, J. Manuf. Sci. Eng. 127 (1) (2005) 23. [40] J. Lee, B.A. Gozen, O.B. Ozdoganlar, Modeling and experimentation of bone drilling forces, J. Biomech. 45 (6) (2012) 1076–1083. [41] R.A. Williams, A study of the drilling process, J. Eng. Ind. 96 (4) (1974) 1207. [42] J. Audy, A study of computer-assisted analysis of effects of drill geometry and surface coating on forces and power in drilling, J. Mater. Process. Technol. 204 (1–3) (2008) 130–138. [43] L.-B. Zhang, L.-J. Wang, X.-Y. Liu, H.-W. Zhao, X. Wang, H.-Y. Luo, Mechanical model for predicting thrust and torque in vibration drilling fibre-reinforced composite materials, Int. J. Mach. Tools Manuf 41 (5) (2001) 641–657. [44] D.E. Brehl, T.A. Dow, Review of vibration-assisted machining, Precis. Eng. 32 (3) (2008) 153–172. [45] Z. Guo, H. Yang, S. Wang, C. Zhou, X. Liu, Adaptive harmony search with bestbased search strategy, Soft. Comput. 22 (4) (2016) 1335–1349. [46] J. Kalivarapu, S. Jain, S. Bag, An improved harmony search algorithm with dynamically varying bandwidth, Eng. Optim. 48 (7) (2015) 1091–1108.