An Iterated Hybrid Local Search Algorithm for Pick-and-Place Sequence Optimization

Abstract

This paper shows the results of our study on the pick-and-place optimization problem. To solve this problem efficiently, an iterated hybrid local search algorithm (IHLS) which combines local search with integer programming is proposed. In the section of local search, the greedy algorithm with distance weight strategy and the convex-hull strategy is developed to determine the pick-and-place sequence; in the section of integer programming, an integer programming model is built to complete the feeder assignment problem. The experimental results show that the IHLS algorithm we proposed has high computational efficiency. Furthermore, compared with the genetic algorithm and the memetic algorithm, the IHLS is less time-consuming and more suitable in solving a large-scale problem.

1. Introduction

The Printed Circuit Board (PCB) is a vital part of electronic products such as smart phones, laptops, and monitors. It contains large amount of components and becomes more and more complicated. In a typical surface mount technology assembly line, PCBs go through a series of processes: solder or adhesive application; component pick-and-place; reflowing; cleaning; testing; and inspection [1]. The pick-and-place process is the most time-consuming procedure. In order to improve production efficiency economically, the pick-and-place sequence optimization is an economical way to reduce production times, and it is the key point to improve production efficiency.

In previous researches, pick-and-place sequence optimization problem has drawn great attentions. There are three sub-problems to be considered in pick-and-place problems: place sequence, feeder assignment and nozzle assignment [2]. Basically, the place sequence (see Section 2.1 steps (2) and (5)) can be considered as the Traveling salesman problem (TSP) or the Vehicle routing problems (VRP) [3,4]. The feeder assignment (see Section 2.1 step (6)) can be viewed as the assignment problem [5]. As to the nozzle assignment, the problem is that different components are assigned to corresponding nozzles, and the object is to minimize the number of nozzle changes. Some researchers do not take the nozzle assignment into consideration in some cases [6,7].

Because of the complexity of the problem, heuristic approaches, such as meta-heuristic approaches and multi-stage heuristics strategy are widely used, while several exact methods can be found in previous research. Meta-heuristic approaches such as genetic algorithm (GA), simulated annealing (SA), particle swarm optimization (PSO) and ant colony optimization (ACO), which play an important role in solving complex problems [8,9,10,11,12], are altered to adapt to solving the problem efficiently [13,14,15,16]. Zhu and Zhang [17] added self-variation behavior to the frog leaping algorithm, which has good accuracy. Multi-stage heuristics strategy, combines the advantages of several algorithm strategies with good performance. Lin and Huang [18] aimed at minimizing total PCB assembly time, and proposed a hybrid genetic algorithm combined with the nearest neighbor search, 2-optimization and genetic algorithm to obtain the optimal component pick-and-place sequence. Similarly, a new mechanism to minimize the assembly time for a multi-headed gantry and high-speed surface mounting technology machine by determining the component assignment to feeder slots was proposed by Han and Seo [19]. The presenting of the heuristic is based on the feeder assignment, which consists of nearest component allocation and globally updated assignment. This algorithm leads to better performance than genetic algorithm or 2-opt swap search. Luo and Liu [20] proposed a two-stage mixed-integer linear programming model and a two-stage problem-solving frame with a hybrid evolutionary algorithm (HEA). The constructive two-stage heuristics is not only determine the set of nozzle types assigned to each head and the total number of assembly cycles, but solve all the sub-problem respectively. As for exact methods, the conventional method is integer programming, and the main measure is modeling and improvement of integer programming model. In research of integer programming models, several researchers focus on the total traveling distance, while others aim at minimizing the assembly completion time. Kumar and Li [21] developed a technique based on integer programming to determine an optimal assignment of feeders as well as an optimal sequence of pickup and placements, furthermore, they obtained optimal solutions based on minimum weight matching and traveling salesman problem. Altinkemer and Kazaz et al. [22] applied an integrated approach which tackles each sub-problems as a single problem simultaneously, and they also presented an $\mathsf{\epsilon }$-approximation algorithm with an $\mathsf{\epsilon }$-error gap for the PCB problem. Ho and Ji [23] formulated several mathematical models for determining the optimal sequence of component placements and the assignment of component types to feeders simultaneously, or for determining the integrated scheduling problem for the sequential pick-and-place machine. The objective of integrated problem is to minimize the total distance traveled by the placement head. Luo and Liu [24] formulated the scheduling optimization problem as a mixed integer linear programming model. And the model was applied to the assembly sequence of placement locations and the assignment of pick-and-place heads for locations so as to minimize the assembly time.

The former researchers widely used effective approaches, such as integer programming, meta-heuristic approaches and hybrid heuristics, to obtain a best solution of pick-and-place sequence optimization problem. These approaches are mainly focused on the optimality of the algorithm results, but paid less attention on time consumption of the algorithm in solving large-scale problems. The aim of present work is to establish an effective multi-stage method with less time-consuming and obtain the optimal solution of the large-scale problem. Considering these restrictions, we propose an integrated hybrid local search (IHLS) algorithm which can meet the time requirements of actual production. Differ from conventional approach which is population based bionic algorithm, the proposed IHLS is a two-stage algorithm which combines the local search method with integer programming method. The main contributions of this paper are: (1) the convex-hull based adjustment strategy is developed to adjust the sequence; (2) a simplified integer programming model of feeder assignment is developed and integrated in the algorithm.

The rest of this paper is organized as follows. In Section 2, the machine structure and PCB assembly process are described, and a simplified mathematical model is proposed. Section 3 is devoted to the description of the IHLS algorithm which combines integer programming with local search. Computational experiments are described in Section 4, which verifies that the algorithm we proposed is superior to genetic algorithm and memetic algorithm in large-scale problems. Finally, conclusions are drawn in Section 5.

2. The Description of Pick-and-Place Sequence Optimization Problem

In this section, the basic structure of pick-and-place machine is described in sufficient detail, and then the PCB assembly processes are introduced. The following description is demonstrating the simplified pick-and-place optimization problem.

2.1. Machine Structure and PCB Assembly Process

As shown in Figure 1 PART I, the main structure of the placement machine is composed of the body frame for fixing other structures, the PCB conveyor for transmitting the PCB, the workbench for fixing and supporting the PCB, the mounting head with a moveable arm for mounting components in x and y directions, and several feeder slots for laying different type of components feeders. As for PCB, as shows in Figure 1 PART III, there is a large amount of points on PCB and each point needs to be mounted with an electrical component. As for feeder, it is a depot for electric components and each feeder only contains one type of electric components.

The assembling process of the PCB by the machine mainly involves the following steps.

(1)

The PCB is transmitted to the workbench by PCB conveyor.

(2)

The PCB is fixed on the workbench, and the coordinates of points on PCB (see Figure 1 PART III) are loaded in the work coordinate system.

(3)

The mounting head moves to the corresponding feeder slots according to a sequence. There are several nozzles on the mounting head, and each nozzle can only grip one component. Thus, the maximum number of components to be carried by the head is the nozzle number (see Figure 1 PART II).

(4)

In order to identify the type of components carried by each nozzle, the mounting head should be scanned by camera before moving to workbench, so as to find the corresponding position on PCB.

(5)

All the components carried by the mounting head are mounted onto PCB.

(6)

The mounting head after mounted moves to the corresponding feeders which are assigned feeder slots, then the head repeat step (3) to step (5) until all the points on PCB are finished. Then, PCB conveyor starts loading the next PCB.

Step (1) and step (2) are the preparing stage in PCB assembly process. Step (3) to step (5) can be defined as a pick-and-place cycle which is the key point of pick-and-place optimization problem. In order to simplify the problem, step (4) is usually ignored in the research. In general, the number of pick-and-place cycles mainly depends on the point number and the nozzle number. The next section will illustrate the simplified problem and notifications.

2.2. The Problem Simplification and Notifications

Step (1) to step (6) are PCB assembly process in real production environment. In researches, the pick-and-place sequence optimization problem is usually simplified while solving the problem. The simplification mainly involves following essential parts.

(1)

Each feeder only contains one component type and can only be assigned to one feeder slot.

(2)

Each point on PCB needs one component and is only allowed to mount once.

(3)

The mounting head needs to pick up components among feeders first, and then the head moves to corresponding points on PCB without scanning procedure by camera.

In addition, the problem assumes that all the nozzles on the head can grip all types of components. Therefore, the exchange of nozzles is not considered in this paper. To better demonstrate the problem, the notifications used in the model are explained as follows.

Parameters:

${\mathit{P}}_{\mathit{X}}={\left(p{x}_1,\dots ,p{x}_i,\dots ,p{x}_N\right)}^T$	The total points’ position in x direction on PCB
${\mathit{P}}_{\mathit{Y}}={\left(p{y}_1,\dots ,p{y}_j,\dots ,p{y}_N\right)}^T$	The total points’ position in y direction on PCB
${\mathit{P}}_{\mathit{T}}={\left(p{t}_1,\dots ,p{t}_i,\dots ,p{t}_N\right)}^T$	The total points need component type on PCB
$N$	The total number of points on PCB
$H$	The number of nozzles
$S$	The total number of feeder slots
$t$	The total number of component types
$n$	The total number of pick-and-place cycles, $n=⎡N/H⎤$
$sx$	The first feeder slot position in x direction
$sy$	The first feeder slot position in y direction
$gap$	The fixed distance between adjacent feeder slots, because all the feeder slots have the same coordinate in y direction, the distance is the difference between two adjacent coordinates in x direction actually
$pD\left(i,j\right)$	The Euclidean distance between point $i$ and point $j$ on PCB
$psD\left(i,p\right)$	The Euclidean distance between point $i$ and feeder slot $p$
$Q_k$	The number of points to be mounted in kth pick-and-place cycle, $Q_k=H$, where $k<n$; $Q_k=N-\left(n-1\right)H$, where $k=n$

The distance $p{D}_{ij}$ between points on PCB can be calculated by formula (1).

$$pD\left(i,j\right)=\sqrt{{\left(p{x}_i-p{x}_j\right)}^2+{\left(p{y}_i-p{y}_j\right)}^2},1\le i\le N,1\le j\le \mathrm{N}$$

(1)

The distance $ps{D}_{ip}$ between point $i$ and feeder slot $p$ can be calculated by formula (2).

$$psD\left(i,p\right)=\sqrt{{\left(p{x}_i-sx-\left(p-1\right)\cdot gap\right)}^2+{\left(p{y}_i-sy\right)}^2},1\le i\le N,1\le 1\le \mathrm{S}$$

(2)

Variables:

$X^k={\left(x_1^k,\dots ,x_i^k,\dots ,x_{Q_k}^k\right)}^T$	$X^k$ is the pick-and-place sequence in kth pick-and-place cycle, $x_i^k\in \left\{1,\dots ,N\right\},1\le i\le Q_k,1\le k\le n$
$a_{gq}\in \left\{0, 1\right\}$	Whether the component type $g$ is assigned to the feeder slot $q$, if $g$ is assigned to slot $q$, $a_{gq}=1$; otherwise, $a_{gq}=0$. $1\le g\le t$, $1\le q\le S$.

2.3. Mathematical Model of the Problem

The aim of this paper is to minimize the total movement distance of mounting head. The object function of the mathematical model of pick-and-place problem is mainly divided into three parts. The first part is component placement distance on PCB, $f_1\left(x\right)$; the second part is the component pick up distance among feeders, $f_2\left(x\right)$; the third part is the movement distance of mounting head between points on PCB and feeders in feeders slots, $f_3\left(x\right)$.

The object function of the model can be described by formulas (3)–(10).

$$F(x)=\min (f_1(x)+f_2(x)+f_3(x))$$

(3)

$$f_1(x)={\displaystyle \sum _{k=1}^n{\displaystyle \sum _{i=2}^{Q_k}pD(x_i^k,x_{i-1}^k)}}$$

(4)

$$f_2(x)={\displaystyle \sum _{k=1}^n{\displaystyle \sum _{i=2}^{Q_k}gap\cdot |{\displaystyle \sum _{q=1}^{S}q}\cdot a_{p{t}_{x_i^k},q}-{\displaystyle \sum _{q=1}^{S}q}\cdot a_{p{t}_{x_{i-1}^k},q}|}}$$

(5)

$$f_3(x)={\displaystyle \sum _{k=1}^{n}psD(x_1^k,{\displaystyle \sum _{q=1}^{S}q\cdot a_{p{t}_{x_i^k},q}})}+{\displaystyle \sum _{k=1}^{n-1}psD(x_H^k,{\displaystyle \sum _{q=1}^{S}q\cdot a_{p{t}_{x_i^{k+1},q}}})}$$

(6)

Subject to

$${\displaystyle \sum _{g=1}^t{a}_{gq}\le 1,q=1,\dots ,S}$$

(7)

$${\displaystyle \sum _{q=1}^S{a}_{gq}=1,g=1,\dots ,t}$$

(8)

$${\displaystyle \sum _{g=1}^t{\displaystyle \sum _{q=1}^S{a}_{gq}}=t}$$

(9)

$${\displaystyle \sum _{k=1}^n{Q}_k=N}$$

(10)

$$a_{gq}\in \{0,1\},g=1,\dots ,t,q=1,\dots ,S$$

(11)

$$x_i^k\in \{1,\dots ,\mathrm{N}\},i=1,\dots ,Q_k$$

(12)

where, formula (7) assures that each feeder slot can be assigned to one type of component feeder at most. Formula (8) assures that each type of component feeder must be assigned to one feeder slot. Formula (9) assures that all the component feeders must be assigned. Formula (10) assures that all the points on PCB must be mounted.

The value of $f_1\left(x\right)$ depends on the place sequence $x_i^k$, and the value of $f_2\left(x\right)$ and $f_3\left(x\right)$ also depend on variable $a_{gq}$. Similarly, $x_i^k$ affects the place sequence $a_{gq}$, which is $a_{p{t}_{x_i^k},q}$. The model of this problem is complicated, which has motivated a lot of researchers to develop kinds of approaches to solve this problem. Generally, it is difficult to solve the problem in an integrated method. Therefore, we propose an algorithm which solves the model in two-stage. The first stage is obtaining the pick-and-place sequence $X^k\left(k=1,\dots ,n\right)$ while the value of $f_1\left(x\right)$ is minimum. The second stage is obtaining the feeder assignment $a_{gq}\left(g=1,\dots ,t,q=1,\dots ,S\right)$ while the value of $f_2\left(x\right)+f_3\left(x\right)$ is minimum.

3. Iterated Hybrid Local Search Algorithm

The iterated hybrid local search algorithm, a two-stage based algorithm, mainly consists of local search and integer programming method. In the first stage, the local search method obtains the pick-and-place sequence $X^k\left(1\le k\le n\right)$ while the value of $f_1\left(x\right)$ is minimum. Based on the obtained $X^k\left(1\le k\le n\right)$, in the second stage, the integer programming aims at obtaining the value of $a_{gq}\left(1\le g\le t,1\le q\le S\right)$ while the value of $f_2\left(x\right)+f_3\left(x\right)$ is minimum. The measures of the local search include greedy strategy with distance weight, adjustment strategy based on convex-hull, termination criteria, and the feeder assignment problem based on the integer programming method. The details of these parts will be explained in the following section. In the algorithm, there are several parameters adopted, and they are explained below.

$f_1\left(x\right)$	The solution of pick-and-place distance in each interaction.
${f^{\prime }}_1\left(x\right)$	The best solution of pick-and-place total distance in iterations.
${X^{\prime }}^k\left(1\le k\le n\right)$	The best solution of pick-and-place sequence in iterations.
${\mathrm{G}}^{\prime }\left(\mathrm{a}\right)$	The best object function value of feeder assignment model in iterations.
$a_{gq}^{\prime }\left(1\le g\le t,1\le q\le S\right)$	The best solution of feeder assignment in iterations.
$\omega \left(i\right)\left(1\le i\le H\right)$	The distance weight which is proposed based on splitting unit circle. $\omega \left(i\right)$ is used to generate initial pick-and-place sequences and it can be obtained by formula (13).

As shown in Figure 2, a circle with a radius of 1 ($r=1$) is divided into H parts, the circumference of circle is $2\pi$, therefore the arc length of each part is $d_{\mathrm{ARC}}=\frac{2\pi }{H}$. Similarly, the angle of each part is $\theta =\frac{2\pi }{H}$, the angle between point i and 1 is $\left(i-1\right)\theta$ which equals to $\frac{2\pi \left(i-1\right)}{H}$. The distance from point $i$ to point 1 is $s_i$ in the circle. The relationship between the distance $s_i$ and $\left(i-1\right)\theta$ can be expressed as an equation, $\cos \left(\left(i-1\right)\theta \right)=\frac{r^2+r^2-s_i^2}{2r}$. The distance $s_i=\sqrt{2-2cos\left(\frac{2\pi \left(i-1\right)}{H}\right)}\left(2\le i\le H\right)$ can be determined. Therefore, $\omega \left(i\right)=\frac{d_{ARC}}{d_{ARC}+s_i}$, the ratio which equals to the distance between the point $i$ and point $i-1$ divide by the sum distance, including the distance between point $i$ and point 1 as well as the distance between point $i$ and point $i-1$, in a unit circle. For example in Figure 2, the angle between point 3 and point 1 is $\frac{2\pi }{3}$, the relationship equation is $\cos \left(\frac{2\pi }{3}\right)=\frac{1+1-s_3^2}{2}$, the distance from point 3 to point 1 is $s_3=\sqrt{2-2\cos \left(\frac{2\pi }{3}\right)}$, while $d_{ARC}=\frac{2\pi }{3}$, then $\omega \left(3\right)=\frac{\pi /3}{\pi /3+s_3}$.

$$\omega \left(i\right)=\frac{\sqrt{2}\pi }{\sqrt{2}\pi +H\sqrt{1-\cos \left(\frac{2\pi \left(i-1\right)}{H}\right)}};2\le i\le H$$

(13)

where, $i$ is the index of nozzles in mounting head.

3.1. The Greedy Strategy with Distance Weight

The pick-and-place sequence initialization is aimed at obtaining initial sequences $X^k\left(k=1,\mathrm{..},n\right)$ which have better qualities. During the iterations of heuristic method, obtaining a better feasible solution contributes to the improvement of efficiency at the beginning. In this section, we propose a greedy strategy with distance weight to initialize the sequences.

After determining $\omega \left(i\right)$, the weighted distance can be represented as formula (14), if point $j$ as the ith sequence in $X^k$.

$${\sigma }_{j,i}^k=\omega \left(i\right)pD\left(j,x_{i-1}^k\right)+\left(1-\omega \left(i\right)\right)pD\left(j,x_1^k\right);2\le i\le Q_k,x_{i-1}^k\in X^k,j\in \left\{1,\dots ,N\right\},j\notin X^k$$

(14)

where ${\sigma }_{j,i}^k$ is the value of weighted distance, $j$ represents the point which is not allocated to $X^k$, and $i$ represents the ith sequence in $X^k$.

$$x_i^k=j, j \mathrm{satisfies} \underset{1\le j\le \mathrm{N},j\notin X^k}{\min }{\sigma }_{j,i}^k$$

(15)

Formula (15) represents that the ith sequence in $X^k$ is the encoding of the point which has a minimum value of ${\sigma }_{j,i}^k$. The pick-and-place sequence initial pseudo code is given below.

Step 1:

Start calculate kth pick-and-place cycle where $k=1$.

Step 2:

If the first iteration involves this process, $x_1^1=1$ and $x_1^k\left(2\le k\le n\right)$ is other random encoding of points; else, $x_1^1=ite{r}_x$ (See Section 3.4) and $x_1^k\left(2\le k\le n\right)$ is other random encoding of points.

Step 3:

In kth pick-and-place cycle, as $x_1^k\left(2\le k\le n\right)$ is confirmed, calculate each ${\sigma }_{j,i}^k$ by formula (14) and let $x_i^k=j$ while the corresponding $j$ satisfies $\underset{1\le j\le N,j\notin X^k}{\min }{\sigma }_{j,i}^k$, from $i=2$ to $i=Q_k$.

Step 4:

If $k=n$, output the initial sequence $X^k\left(1\le k\le n\right)$; otherwise, $k=k+1$ and repeat from step 3 to step 4.

For example, the total number of points on PCB is 6, the coordinates of points are shown in

Table 1. The example of greedy strategy with distance weight.

No.	$\mathit{x}$	$\mathit{y}$
1	16	2
$2$	18	2
$3$	19	12
$4$	17	14
5	11	8
6	15	4

. The number of nozzles is 6, the pick-and-place sequence cycle number is 1, as shown in Figure 2. In order to initialize $X^1$, $\omega \left(i\right)\left(2\le i\le 6\right)$ can be determined, which are $\omega \left(2\right)=0.5115$, $\omega \left(3\right)=0.3767$, $\omega \left(4\right)=0.3435$, $\omega \left(5\right)=0.3765$, $\omega \left(6\right)=0.5108$. As $x_1^1=1$, the next sequence $x_2^1$ need to be determined based on the values of ${\sigma }_{2,2}^1$, ${\sigma }_{3,2}^1$, ${\sigma }_{4,2}^1$, ${\sigma }_{5,2}^1$ and ${\sigma }_{6,2}^1$. The values of above items are ${\sigma }_{2,2}^1=2.0$, ${\sigma }_{3,2}^1=10.4$, ${\sigma }_{4,2}^1=12.0$, ${\sigma }_{5,2}^1=7.8$ and ${\sigma }_{6,2}^1=2.2$. According to formula (15), $x_2^1=2$, which is the minimun among ${\sigma }_{2,2}^1,{\sigma }_{3,2}^1,{\sigma }_{4,2}^1,{\sigma }_{5,2}^1\mathrm{and}{\sigma }_{6,2}^1$. The next squence $x_3^1$ depends on the value of ${\sigma }_{3,2}^1$, ${\sigma }_{4,2}^1$, ${\sigma }_{5,2}^1$ and ${\sigma }_{6,2}^1$. In this way, the sequence of $X^1$ can be determined. This strategy is aimed at initializing the sequence $X^k$ so that it is a better feasible solution, but the length of each sequence $X^k$ is not the minimum. So the adjusting strategy is needed to get a sequence of minimum distance.

3.2. The Adjusting Strategy Based on Convex-Hull

To adjust the sequence, each initial pick-and-place sequence can be regarded as a small size traveling salesman problem (TSP). For the traveling salesman problem, there are several classifications and extensions, such as symmetric traveling salesman problem (STSP), asymmetric traveling salesman problem (ATSP) and multiple traveling salesman problem (MTSP). Laszlo Barna Iantovics [25] proposed a novel universal metric and solved STSP by intelligent cooperative multi-agent systems. Furqan Hussain Essani [26] presented an algorithm combined colored petri nets and guaranteed a feasible solution to the MTSP with asymmetric cost matrix. Differ from above methods, the TSP in this section needs several simple measures to guarantee the efficiency of the algorithm and obtain a better feasible solution. Based on convex-hull algorithms [27], an approximate method is developed, and the details are explained below.

$$\rho \left(a,b,c\right)=pD\left(a,b\right)+pD\left(a,c\right)-pD\left(b,c\right),a,b,c\in X^k$$

(16)

where $\rho \left(a,b,c\right)$ represents the additional length, if line $bc$ is replaced by line $ab$ and line $ac$.

The pick-and-place sequence adjustment pseudo code of each cycle is given below.

Step 1:

Let $A={\left(x_1^k,x_2^k,x_3^k\right)}^{\mathrm{T}}$, $U={\left(x_4^k,\dots ,x_{Q_k}^k\right)}^{\mathrm{T}}$.

Step 2:

Select each item $x_p^k$ from $U$, $x_p^k\in U$; Select each pair of items $x_p^k,x_{p+1}^k$ from $A$, $x_p^k,x_{p+1}^k\in A$. Calculate the value of $\rho \left(x_p^k,x_q^k,x_{q+1}^k\right)$.

Step 3:

Find the corresponding $x_p^k$, where $x_p^k$ satisfies $\underset{x_p^k\in U and x_q^k,x_{q+1}^{\mathrm{k}}\in A}{\min }\rho \left(x_p^k,x_q^k,x_{q+1}^k\right)$; put $x_p^{\mathrm{k}}$ into $A$ between $x_q^k$ and $x_{q+1}^k$, remove $x_p^k$ from $U$.

Step 4:

If $U=\varphi$, turn to step 7; otherwise, repeat from step 2 to step 3.

Step 5:

Update $X^k$, let $X^k=A$.

For example, as shown in Figure 3, $A=\left(1,2,6\right)$, and $U=\left(3,4,5\right)$. The dashed lines represent the additional length when the point is added to $A$, the solid lines represent the length between two adjacent points in $A$ sequence. Each additional length of item in $U$ ($x_p=5,4\mathrm{or}3$) and line 2–6 in $A$ ($x_q=2,x_{q+1}=6$) can be obtained, which are $\rho \left(5,6,2\right)$, $\rho \left(3,6,2\right)$ and $\rho \left(4,6,2\right)$. The smallest value of $\rho$ is $\rho \left(3,6,2\right)$, and then $A=\left(1,2,3,6\right)$, $U=\left(4,5\right)$.

During the sequence adjustment, the additional length of each adjusting measure (step 3) is considered, and only the adjustment of minimum length will be accepted in this strategy. In other words, all cases of adjustment can be accepted based on formula (16). Therefore, the length of initial $X^k$ is minimum after the adjusting strategy.

3.3. Feeder Slots Selection

The feeder slots selection is based on the pick-and-place sequence $X^k\left(1\le k\le n\right)$. In order to improve efficiency, a new feeder assignment integer programming model is built in this section, and the model is solved by integer programming method. This model is modified based on the mathematical model (Section 2), replaced several variables, and retained the original constraints. As $X^k\left(1\le k\le n\right)$ is determined, we generate parameters $lin{k}_{fg}$ and $boun{d}_r$.

$lin{k}_{fg}$	The number of links between component type $f$ and $g$, where $1\le f,g\le t$ and $f\ne g$.
$boun{d}_r$	The first or last point in each pick-and-place cycle $X^k$, where $1\le r\le 2n$.

Because $x_i^k\in X^k$ is the encoding of points, $p{t}_{x_i^k}$ represents the corresponding component type of $x_i^k$. As $X^k=\left(x_1^k,x_2^k,\dots ,x_i^k,\dots ,x_{Q_k}^k\right)$ is known, the corresponding component types of $X^k$ are $\left(p{t}_{x_1^k},p{t}_{x_2^k},\dots ,p{t}_{x_i^k},\dots ,p{t}_{x_{Q_k}^k}\right)$, which $p{t}_{x_i^k}\in \left\{1,\dots ,t\right\}$. In the component type vector $\left(p{t}_{x_1^k}, p{t}_{x_2^k},\dots ,p{t}_{x_i^k},\dots ,p{t}_{x_{Q_k}^k}\right)$, $p{t}_{x_1^k},p{t}_{x_2^k}$ is a link between component type $p{t}_{x_1^k}$ and $p{t}_{x_2^k}$. The parameter $lin{k}_{fg}$ is the total number of links between component type $f$ and $g$, and it can be obtained by result $X^k$. Initialize $lin{k}_{fg}=0 \left(1\le f,g\le t\right)$, then $lin{k}_{fg}$ can be obtained from $\left(p{t}_{x_1^k},p{t}_{x_2^k},\dots ,p{t}_{x_i^k},\dots ,p{t}_{x_{Q_k}^k}\right)$ as formula (17).

$$lin{k}_{p{t}_{x_i^k},p{t}_{x_{i+1}^k}}=lin{k}_{p{t}_{x_i^k},p{t}_{x_{i+1}^k}}+1,p{t}_{x_i^k}\ne p{t}_{x_{i+1}^k},1\le i\le Q_k-1,1\le k\le n$$

(17)

For example, the component type vector is $\left(p{t}_2,p{t}_1,p{t}_4,p{t}_4,p{t}_2,p{t}_2\right)$, the corresponding links can be obtained as following

Table 2. The corresponding $lin{k}_{fg}$.

$\mathit{l}\mathit{i}\mathit{n}{\mathit{k}}_{\mathit{f}\mathit{g}}$	$\mathit{p}{\mathit{t}}_1$	$\mathit{p}{\mathit{t}}_2$	$\mathit{p}{\mathit{t}}_3$	$\mathit{p}{\mathit{t}}_4$
$p{t}_1$	0	0	0	1
$p{t}_2$	1	0	0	0
$p{t}_3$	0	0	0	0
$p{t}_4$	1	1	0	0

.

As for $boun{d}_r$, formula (18) demonstrates the relationship between $boun{d}_r$ and $X^k$.

$$\{\begin{array}{c}boun{d}_{2k-1}=x_1^k, 1\le k\le n\\ boun{d}_{2k}=x_{Q_k}^k, 1\le k\le n\end{array}$$

(18)

In the second stage, the object function should be $f_2\left(x\right)+f_3\left(x\right)$ (formulas (5) and (6) in Section 2.1). It can be expressed as ${\displaystyle \sum }_{k=1}^n{\displaystyle \sum }_{i=1}^{Q_k}gap\left|{\displaystyle \sum }_{q=1}^{S}q·a_{p{t}_{x_i^k,q}}-{\displaystyle \sum }_{q=1}^{S}q·a_{p{t}_{x_{i-1}^k,q}}\right|+{\displaystyle \sum }_{k=1}^{n-1}psD\left(x_H^k,{\displaystyle \sum }_{q=1}^{S}q·a_{p{t}_{x_i^{k-1},q}}\right)$. ${\displaystyle \sum }_{q=1}^{S}q·a_{p{t}_{x_i^k,q}}$ represents the corresponding component type, when $x_i^k$ is obtained in first stage. Based on above works, ${\displaystyle \sum }_{\mathrm{q}=1}^{\mathrm{S}}q·a_{fg}·lin{k}_{fg}$ replaces ${\displaystyle \sum }_{q=1}^{S}q·a_{p{t}_{x_i^k,q}}$, and $boun{d}_r$ replaces $x_H^k$. Without modifying of constraints, the new local integer programming model is described as follows:

$$G(a)=\min {\displaystyle \sum _{f=1}^t{\displaystyle \sum _{g=1}^{t}gap\cdot lin{k}_{fg}\left|{\displaystyle \sum _{q=1}^{S}q{a}_{fq}}-{\displaystyle \sum _{q=1}^{S}q{a}_{gq}}\right|}}+{\displaystyle \sum _{r=1}^{2n}{\displaystyle \sum _{q=1}^S{\displaystyle \sum _{g=1}^t{a}_{gq}}\cdot psD(boun{d}_r,q)}}$$

(19)

Subject to

$${\displaystyle \sum _{g=1}^t{a}_{gq}\le 1,g=1,\dots ,t}$$

$${\displaystyle \sum _{q=1}^t{a}_{gq}=1,g=1,\dots ,t}$$

$${\displaystyle \sum _{g=1}^t{\displaystyle \sum _{q=1}^S{a}_{gq}}=t}$$

$$a_{gq}\in \{0,1\};g=1,\dots ,t,q=1,\dots ,S$$

Binary variable $a_{gq}$ represents whether the ith component type feeder is in the slot $q$, and ${\displaystyle \sum }_{q=1}^{S}q\cdot a_{gq}$ represents the corresponding sequence of $g$th component type. Therefore, $\left|{\displaystyle \sum }_{q=1}^{S}q\cdot a_{fq}-{\displaystyle \sum }_{q=1}^{S}q\cdot a_{gq}\right|$ indicates the number of feeders between slot $f$ and $g$. While ${\displaystyle \sum }_{f=1}^t{\displaystyle \sum }_{g=1}^{t}gap·lin{k}_{fg}\left|{\displaystyle \sum }_{q=1}^{S}q{a}_{fq}-{\displaystyle \sum }_{q=1}^{S}q{a}_{gq}\right|$ demonstrates the total distance of the mounting head movements among feeders, it is the same as $f_2\left(x\right)$ in Section 2.3. $boun{d}_r\left(1\le r\le 2n\right)$ is encoded of first or last point in each sequence, while $psD\left(boun{d}_r,q\right)$ is the distance between the point $boun{d}_r$ and feeder slot $q$. Therefore, $psD\left(boun{d}_r,q\right)·{\displaystyle \sum }_{g=1}^t{a}_{gq}$ represents the distance between point $boun{d}_r$ and the corresponding feeder type slot. ${\displaystyle \sum }_{r=1}^{2n}{\displaystyle \sum }_{q=1}^S{\displaystyle \sum }_{g=1}^{t}psD\left(boun{d}_r,q\right)·{\displaystyle \sum }_{g=1}^t{a}_{gq}$ demonstrates the total moving distance of mounting head between points on PCB and feeders in slots, which is the same as $f_3\left(x\right)$ in Section 2.3.

The object function is to minimize the total distance of mounting head based on the pick-and-place sequence $X^k\left(1\le k\le n\right)$. And the constraints of this model are as same as the model in Section 2.3.

3.4. Initialization Update and Termination Criteria

In this algorithm, we obtain a feasible solution including pick-and-place sequence $X^k\left(1\le k\le n\right)$ and feeder assignment $a_{gq}\left(g=1,\dots ,t,q=1,\dots ,S\right)$. To improve the quality of solution, selecting $X^k$ based on the maximum length and resetting sequence $X^k$ are crucial iteration strategy. During each iteration, $X^k$ needs to be recorded if the value of the length is maximum, and this sequence needs to be reset first in next iteration, which is $x_1^1=ite{r}_x$ (see Section 3.1). $ite{r}_x$ can be generated by formula (20).

$$ite{r}_x=x_1^k, k satisfies \underset{1\le k\le n }{\max }{\displaystyle \sum }_{i=2}^{Q_k}pD\left(x_i^k,x_{i-1}^k\right)$$

(20)

The IHLS algorithm is terminated when the number of iterations equals to the maximum number of iterations. The flowchart of the IHLS (Figure 4) and pseudo code are given below.

Step 1:

Initialize the parameters.

Step 2:

Initialize the pick-and-place sequence $X^k$ with $\omega \left(i\right)$ according to greedy strategy, where $1\le k\le n$.

Step 3:

Adjust the pick-and-place sequence $X^k$ by adjustment strategy, where $1\le k\le n$.

Step 4:

Calculate the fitness value $f_1\left(x\right)$.

Step 5:

If the best solution $f_1^{\prime }\left(x\right)>f\left(x\right)$, turn to step 6; otherwise, repeat from step 3 to step 5.

Step 6:

Update the record best solution $f_1^{\prime }\left(x\right)$, let $f_1^{\prime }\left(x\right)=f\left(x\right)$ and ${X^{\prime }}^k=X^k$.

Step 7:

According to ${X^{\prime }}^k$, relevant parameters ($lin{k}_{fg}$ and $boun{d}_r$) are prepared and input to feeder assignment model.

Step 8:

Update the best solutions $G^{\prime }\left(a\right)$ and $a_{gq}^{\prime }$, let $G^{\prime }\left(a\right)=G\left(a\right)$ and $a_{gq}^{\prime }=a_{gq}$.

Step 9:

If termination criteria was satisfied, turn to step 10; otherwise, repeat from step 3 to step 8.

Step 10:

Output the best solutions.

4. Experimental Results and Analysis

In order to verify the effectiveness of IHLS algorithm, three other test algorithms are selected for comparison. The three baseline algorithms are genetic algorithm (GA), memetic algorithm and integer programming method. The algorithms are implemented in MATLAB and run with the 2.8 GHz computer with 16 GB RAM. The integer programming models are solved by Yalmip with Cplex1263.

4.1. Comparison with Exact Algorithm

As for integer programming, because integer programming method can only solve very small-scale problems, the place sequence integer programming model and several small instances except instance 4 are generated based on the research of Ho and Ji [2].

The IHLS is executed on very small-scale instances and the results are obtained by integer programming method, Cplex. Three small-scale instances with 6 feeder slots are given in

Table 3. Comparison of the IHLS results with the optimal solutions.

(PN, TN)	No.	Object Function Value (mm)			CUP Times
		${\mathit{D}}_{\mathit{I}\mathit{H}\mathit{L}\mathit{S}}$	${\mathit{D}}_{\mathit{C}\mathit{p}\mathit{l}\mathit{e}\mathit{x}}$	$\frac{{\mathit{D}}_{\mathit{I}\mathit{H}\mathit{L}\mathit{S}}-{\mathit{D}}_{\mathit{C}\mathit{p}\mathit{l}\mathit{e}\mathit{x}}}{{\mathit{D}}_{\mathit{C}\mathit{p}\mathit{l}\mathit{e}\mathit{x}}}\times 100\%$	$\mathit{C}{\mathit{T}}_{\mathit{I}\mathit{H}\mathit{L}\mathit{S}}$	$\mathit{C}{\mathit{T}}_{\mathit{c}\mathit{p}\mathit{l}\mathit{e}\mathit{x}}$
(6, 5)	1	280.12	263.37	6.45	0.35	0.27
(8, 6)	2	409.41	383.71	6.78	0.98	0.83
(10, 6)	3	598.34	566.02	5.65	2.21	2.09
(82, 3)	4	3766	Unsolved	-	22.3	-

. The results show that IHLS is near to the optima as the average gap is 6.29%. While the instance 4 (the nozzle number of instance 4 is 1) shows that the integer programming method cannot solve the problem with a larger scale.

4.2. Comparison of the Heuristics in Different Size Instances and Analysis

In this experiment, the performance of the algorithm is tested on 9 PCB instances. These 9 PCB instances are generated based on industrial data.

(1) Set the model parameters

The model parameters, such as $N$, $t$, $n$, ${\mathit{P}}_{\mathit{X}}$, ${\mathit{P}}_{\mathit{Y}}$ and ${\mathit{P}}_{\mathit{Y}}$, are based on the test instances, and other parameters, such as $H$, $S$, $\mathrm{s}x$, $\mathrm{s}y$, $gap$ and $Q_k$, are based the pick-and-place machine parameters. The parameters are introduced in

Table 4. The size of instances.

No.	N	t	H	n	S	(sx, sy)	$\mathbf{gap}$	${\mathbf{Q}}_{\mathit{k}}$
1	82	3	10	9	60	(−120, −72)	12	(10, 10, …, 3)
2	100	4	10	10	60	(−120, −72)	12	(10, 10, …, 10)
3	145	5	10	15	60	(−120, −72)	12	(10, 10, …, 5)
4	163	4	10	17	60	(−120, −72)	12	(10, 10, …, 3)
5	184	5	10	19	60	(−120, −72)	12	(10, 10, …, 4)
6	200	5	10	20	60	(−120, −72)	12	(10, 10, …, 10)
7	240	6	10	24	60	(−120, −72)	12	(10, 10, …, 10)
8	264	7	10	27	60	(−120, −72)	12	(10, 10, …, 4)
9	300	9	10	30	60	(−120, −72)	12	(10, 10, …, 10)

.

The total number of components ($N$) ranges from 82 to 300 and the component type (t) ranges from 3 to 9. The solving complexity of instances mainly depend on the number of components. With the increasing of $N$, the difference between algorithms performance can be more and more significant.

(2) Comparison with the heuristics

In the GA algorithm, the chromosome structure is decimal code, and each chromosome includes PCB placement sequence part and feeder assignment part. The chromosome structure can be found in Kulak, Yilmaz and Günther [13]. The population size is 2000, the crossover rate is 90%, and the mutation rate is 10%. As for memetic algorithm, based on Carlos’ research [28], the memetic algorithm used GA frame in population-based search, and the greedy strategy is used in local search part. The parameters of the memetic algorithm are the same as GA’s. The three algorithms, IHLS, GA and memetic algorithm, solve these instances based on the mathematical model in Section 2. The optimization results of each instance are the average values of 5 independent runs. The solutions are presented in

Table 5. Experimental results of different size data.

No.			GA		IHLS		Memetic Algorithm
	${\mathbf{D}}_{\mathbf{best}}$	$\mathbf{Minimum}\mathbf{CT}$	${\mathbf{D}}_{\mathbf{GA}}\left(\mathbf{mm}\right)$	$\mathbf{C}{\mathbf{T}}_{\mathbf{GA}}$	${\mathbf{D}}_{\mathbf{IHLS} }\left(\mathbf{mm}\right)$	$\mathbf{C}{\mathbf{T}}_{\mathbf{IHLS}}$	${\mathbf{D}}_{\mathbf{MA}}\left(\mathbf{mm}\right)$	$\mathbf{C}{\mathbf{T}}_{\mathbf{MA}}$
1	3511	22.3	3542	1534.0	3766	22.3	3511	1485.7
2	5509	29.5	5377	1512.4	5844	29.5	5509	1508.2
3	7473	52.4	7655	1743.5	7740	52.4	7473	1485.0
4	8166	65.0	8352	1744.0	8479	65.0	8166	1846.8
5	8813	78.4	10,439	1766.8	8813	78.4	9430	2020.5
6	9882	93.4	11,541	1923.4	9882	93.4	10,249	2195.0
7	11,218	124.2	14,745	2043.8	11,218	124.2	11,725	2596.0
8	13,591	174.2	17,610	2142.0	13,591	174.2	13,837	2802.0
9	14,854	247.0	20,745	2689.4	14,854	247.0	15,007	3310.7

. ${\mathrm{D}}_{\mathrm{GA}}$, ${\mathrm{D}}_{\mathrm{IHLS}}$ and ${\mathrm{D}}_{\mathrm{MA}}$ represent the total distance of GA, IHLS and memetic algorithm respectively. ${\mathrm{CT}}_{\mathrm{GA}}$, ${\mathrm{CT}}_{\mathrm{IHLS}}$ and ${\mathrm{CT}}_{\mathrm{MA}}$ represent the CPU times of each algorithm.

(3) Experimental results analysis

There are different iteration strategies in three algorithms. The three algorithms obtain the optimal solution when the iteration achieves the maximum iteration number. But the solutions are different among the three algorithms. As shown in the second, forth, sixth and eighth column in Table 5, from instance 1 to instance 4 except instance 2, the minimum distance value obtained by memetic algorithm can be found. As for GA and IHLS, the two algorithms’ value of distance is larger than memetic algorithm’s. From instance 5 to instance 9, the distance of IHLS is the minimum, and memetic algorithm’s value of distance is close to that of IHLS algorithm. However, as for GA, the distance is much larger than IHLS’s and memetic algorithm’s.

Considering the CPU times of each algorithm in Figure 5, IHLS is the most efficient one among the three algorithms. As shown in Figure 5, the swarm based algorithms, GA and memetic algorithm, are more time-consuming. With the increasing of instance size, the CPU time of memetic algorithm is getting larger, and memetic algorithm become the most time-consuming one.

In respect of optimality, IHLS algorithm has an advantage in large-scale instances, while the solution of IHLS is slightly inadequate in small one. In respect of CPU times, IHLS algorithm has a great advantage in both small-scale and large-scale instances. To show the convergence of different heuristics, the stopping criteria is set as 2000 iterations. The convergence plot comparisons of the three algorithms are shown in Figure 6 and Figure 7. Figure 6 shows the convergence plot in solving instance 1, 2 and 3 which are small-size data, while Figure 7 shows the convergence plot in solving instance 7, 8 and 9 which are large-size data. More analysis details will be explained in the following.

$$\tau =\frac{D_{HILS}-D_{best}}{D_{best}}\times 100\%$$

(21)

To measure IHLS performance, a relative quality measure $\frac{D_{HILS}-D_{best}}{D_{best}}\times 100\%$, where $D_{best}$ is the best among the solutions obtained by the three algorithms and $D_{HILS}$ is IHLS. The closer to 0 the measure value $\tau$ is, the better the according algorithm will be.

As shown in Figure 6, from instance 1 to instance 3, the number of points $N$ < 150, there is a flaw in IHLS’ solutions. GA obtained minimum object function value which is the best feasible solution in instance 2, memetic algorithm obtained the best feasible solution in instance 1 and instance 3. In Figure 6 (instance 1), the best solution is 3511 mm, the feasible solution of IHLS is 3766 mm, $\tau =7.26\%$. Similarly, the measure value $\tau$ is 8.68% in Figure 6 (instance 2). While the measure value $\tau =3.57\%$ in Figure 6 (instance 3), which becomes smaller. Although IHLS lacks good performance in Figure 6, the difference between the IHLS’ object function and the best solution is less than 700. As shown in Figure 7, from instance 7 to instance 9, the number of points $N>200$, IHLS’s solution is the best among three heuristics. Obviously, the measure value $\tau$ is $0$ in Figure 7. In Figure 7, memetic algorithm also has a good performance, but the quality of solution is slightly worse than that of IHLS. The difference between memetic algorithm’s and IHLS’ solution values are 507, 240 and 157 in instance 7, 8 and 9 respectively, while the difference of GA are more than 3000. Figure 6 and Figure 7 show that IHLS lacks good performance in solving small-scale problem, while with the increasing of points number, it demonstrates that IHLS has a better performance in solving large-scale problem.

From instance 1 to 9, both IHLS and memetic algorithm achieve a better convergence value in the first iteration, while the GA’s first convergence value is worse. The solutions of IHLS and memetic algorithm converges relatively rapidly. In addition, IHLS is more suitable for real-time application in operation workshop because of its higher computational efficiency. Considering the time-consuming and optimality, IHLS is superior to memetic algorithms.

The results of instance 9 are demonstrated in

Table 6. Results of IHLS in instance 9.

$\mathit{k}$	${\mathbf{X}}^{\mathbf{k}}$	${\mathit{a}}_{\mathit{g}\mathit{q}}$	Corresponding Feeder Type	Corresponding Feeder Coordinates
1	$X^1=\left(146, 219, 13, 27, 185, 224, 284, 259, 215, 233\right)$	$a_{1,27}=1$	100004652	(192, −72)
2	$X^2=\left(93, 72, 113, 165, 177, 204, 87, 31, 89, 152\right)$	$a_{2,23}=1$	100004697	(144, −72)
3	$X^3=\left(108, 237, 153, 174, 101, 69, 210, 140, 126, 62\right)$	$a_{3,24}=1$	100004795	(156, −72)
4	$X^4=\left(105, 48, 58, 186, 143, 129, 77, 57, 231, 176\right)$	$a_{4,26}=1$	100004883	(180, −72)
5	$X^5=\left(66, 104, 128, 268, 214, 184, 276, 280, 239, 238\right)$	$a_{5,16}=1$	100005022	(16, −72)
6	$X^6=\left(63, 234, 179, 131, 65, 75, 183, 132, 47, 264\right)$	$a_{6,3}=1$	100006424	(−96, −72)
7	$X^7=\left(144, 74, 59, 213, 158, 159, 98, 99, 270, 45\right)$	$a_{7,22}=1$	100008025	(132, −72)
8	$X^8=\left(245, 29, 226, 187, 241,251, 190, 188, 253, 260\right)$	$a_{8,25}=1$	100030843	(168, −72)
9	$X^9=\left(244, 271, 265, 200, 203, 266, 209, 212, 211, 230\right)$	$a_{8,15}=1$	100031371	(48, −72)
10	$X^{10}=\left(293, 235, 269, 205, 206, 267, 279, 202,236, 199\right)$
11	$X^{11}=\left(141, 156, 68, 21, 207, 180, 168, 110, 246, 111\right)$
12	$X^{12}=\left(232, 181, 295, 257, 258, 229, 208, 248, 17, 262\right)$
13	$X^{13}=\left(26, 247, 220, 194, 242, 243, 193, 217, 256, 10\right)$
14	$X^{14}=\left(223, 197, 282, 281, 196, 218, 263, 225, 164, 161\right)$
15	$X^{15}=\left(44, 95, 52, 119, 147, 162, 198, 27, 49, 71\right)$
16	$X^{16}=\left(22, 7, 189, 123, 122, 83, 9, 216, 135, 138\right)$
17	$X^{17}=\left(252, 11, 255, 191, 227, 250, 249, 221, 182, 254\right)$
18	$X^{18}=\left(195, 137, 134, 86, 170, 125, 80, 81, 192, 8\right)$
19	$X^{19}=\left(240, 32, 222, 149, 155, 107, 36, 228, 167, 16\right)$
20	$X^{20}=\left(173, 92, 12, 201, 171, 150, 116, 54, 24, 261\right)$
21	$X^{21}=\left(291, 286, 33, 61, 6, 172, 117, 103, 35, 145\right)$
22	$X^{22}=\left(287, 274, 37, 67, 28, 136, 114, 91, 120, 154\right)$
23	$X^{23}=\left(4, 70, 55, 163, 2, 112, 30, 157, 300, 299\right)$
24	$X^{24}=\left(42, 73, 96, 151, 40, 94, 3, 133, 277, 273\right)$
25	$X^{25}=\left(46, 76, 51, 148, 34, 64, 5, 169, 296, 292\right)$
26	$X^{26}=\left(19, 88, 50, 127, 90, 115, 43, 121, 290, 283\right)$
27	$X^{27}=\left(14, 288, 139, 23, 100, 18, 142, 97, 102, 294\right)$
28	$X^{28}=\left(285, 289, 78, 85, 20, 124, 53, 106, 15, 130\right)$
29	$X^{29}=\left(297, 278, 84, 79, 60, 160, 56, 118, 41, 175\right)$
30	$X^{30}=\left(298, 38, 82, 166, 1, 39, 178, 109, 25, 275\right)$

. In the table, the second column shows the pick-and-place sequence in each cycle, and other columns show the assignment of feeders. The path of instance 9 is shown in Figure 8.

5. Conclusions

In this paper, we deal with the pick-and-place optimization problem efficiently by a new algorithm (IHLS) combining local search with integer programming. In the section of local search, the greedy with distance weight strategy and the convex-hull strategy are developed to determine the pick-and-place sequence; in the section of integer programming, an integer programming model is built to solve the feeder assignment. The experimental result shows that the IHLS we proposed has a higher computational efficiency. The IHLS is less time-consuming and more suitable in solving large scale problems. In our future work, in order to improve the accuracy of IHLS in small size data, some other solving techniques could be integrated and the mathematical model should be improved.