Application of Harmony Search to Vehicle Routing

: A phenomenon-inspired meta-heuristic algorithm, harmony search, imitating music improvisation process, is introduced and applied to vehicle routing problem, then compared with one of the popular evolutionary algorithms, genetic algorithm. The harmony search algorithm conceptualized a group of musicians together trying to search for better state of harmony. This algorithm was applied to a test traffic network composed of one bus depot, one school and ten bus stops with demand by commuting students. This school bus routing example is a multi-objective problem to minimize both the number of operating buses and the total travel time of all buses while satisfying bus capacity and time window constraints. Harmony search could find good solution within the reasonable amount of time and computation.


INTRODUCTION
Traffic researchers and engineers sometimes cope with situations where optimization is needed to find the optimal solution out of huge amount of combinatorial solutions. Traditionally efficient enumeration methods such as branch and bound (B&B) technique have been used for supporting these optimal decisions. However, their computational disadvantages such as requiring huge amount of computation and memory made people rely on another type of methodology, that is, evolutionary or meta-heuristic algorithms.
The basic ideas of existing evolutionary and metaheuristic algorithms are motivated by natural phenomena. For example, the evolutionary algorithms [1][2][3] and the genetic algorithm (GA) [4,5] are inspired by biological evolutionary process; tabu search [6] and ant algorithm [7] from animal's behavior; and simulated annealing [8] from physical annealing process.
Harmony Search (HS) algorithm has been recently developed in an analogy with music improvisation process where musicians in an ensemble continue to polish their pitches in order to obtain better harmony [9] . The HS algorithm has been successfully applied to various real-world combinatorial optimization problems such as truss structure design, pipe network design, pump switching, and hydrologic parameter calibration [10][11][12][13] . Especially, for the pump switching problem which has 2 40 ( 1.1× 10 12 ) different combinations, HS found the less energy consuming operation than GA or even than B&B methods which were implemented as an IBM subroutine or MS-Excel Solver. Consequently, the HS algorithm provides a possibility of success in a combinatorial optimization problem in traffic engineering field.
In order to demonstrate the searching ability of HS, it is applied to a school bus routing problem (SBRP) which is a practical optimization problem. From a school's perspective, the SBRP aims to provide students with an efficient and equitable transportation service [14] . Many approaches have been developed to solve the SBRP [14][15][16] . This problem falls into a large class of problems called vehicle routing problem (VRP), in which a set of vehicles provides pickup, delivery or simply a service to customers dispersed in an area.
Generally, the wider the problem scope is, the harder the problem becomes in term of the solution technique such as B&B method. Because an exact method may take a very long time to obtain the optimal solution, people have turned to evolutionary or metaheuristic algorithms such as GA that do not necessarily find the optimal solution, but tend to find good solutions within a reasonable amount of time [17] . Pattnaik et al. [18] proposed an optimization model of minimizing overall cost (both the operator's cost and user's) while determining a route configuration with a set of urban bus transit routes and associated frequencies using GA. In their approach, a set of candidate routes is first generated, and then GA is employed to find the optimal one. Chien et al. [19] introduced a GA model to the optimization of bus route and the corresponding headway while minimizing total cost, subject to geography, capacity, and budget constraints. Their results were validated by comparing with those obtained from exhaustive search algorithms. Jung and Haghani [20] used GA to solve a multi-vehicle pickup and delivery problem with time window. They formulated the problem as a mixed-integer linear program, with an objective function to minimize the total cost, consisting of the fixed cost of the vehicles, routing cost, and customer inconvenient cost. In this study, a recently-developed music phenomenoninspired HS algorithm is applied to the test network of SBRP for the first time and compared with one of the most popular algorithms, GA.

HARMONY SEARCH ALGORITHM
From the idea that existing evolutionary or metaheuristic algorithms imitate natural, scientific, or behavioral phenomena around us, a new algorithm can be conceptualized from a music performance process (say, a jazz improvisation) involving searching for a better harmony. Just like music improvisation seeks a best state (fantastic harmony) determined by aesthetic estimation, optimization process seeks a best state (global optimum) determined by objective function evaluation; Just like aesthetic estimation is determined by the set of the pitches played by ensembled instruments, function evaluation is determined by the set of the values assigned for decision variables; Just like aesthetic sound quality can be improved practice after practice, objective function value can be improved iteration by iteration. The analogy between improvisation and optimization is shown in Fig. 1 x = {500, 600, 700}). If the saxophonist toots the note Do, the double bassist plucks Mi, and the guitarist plucks Sol, their notes together make a new harmony (Do, Mi, Sol). If this new harmony is better than existing harmony, the new harmony is kept. Likewise, the new solution vector (100mm, 300mm, 500mm) generated in optimization process is kept if it is better than existing harmony in terms of objective function value. Just as the harmony quality is enhanced practice after practice, the solution quality is enhanced iteration by iteration.
According to the above algorithm concept, the HS algorithm consists of the following five steps: parameter initialization; harmony memory initialization; new harmony improvisation; harmony memory update; and termination criterion check.
Parameter initialization: In the first step, the optimization problem is specified as follows: . However, there is no sorting required in this vehicle routing problem because each variable value represents just node number rather than physical amount such as node demand); N is the number of decision variables (= number of music instruments); and K is the number of candidate values for the discrete decision variables (= number of pitches for each instrument).
The HS algorithm parameters are also specified in this step: HMS (harmony memory size; = number of simultaneous solution vectors in harmony memory), HMS (harmony memory considering rate), PAR (pitch adjusting rate), and number of improvisations (= number of objective function evaluations). These algorithm parameters are explained in the following steps.

Harmony memory initialization: In
Step 2, the Harmony Memory (HM) as shown in Equation 3, is crammed with as many randomly generated solution vectors as the size of the HM (i.e., HMS). (3)
Random selection: As a musician plays any pitch within the instrument range (for example, {Do, Re, Mi, Fa, Sol, La, Si} in Fig. 2), the value of decision variable i x′ is randomly chosen within the value range i X . HM Consideration. As a musician plays any pitch out of the preferred pitches in his/her memory (for example, {Do, Mi, Do, Sol, Do} in Fig. 3), the value of decision variable i x′ is chosen from any pitches stored ) with a probability of HMCR (0 HMCR 1) while it is randomly chosen with a probability of (1-HMCR) in random selection process as previously described.
Pitch Adjustment. Once one pitch is obtained in HM consideration, a musician can further adjust the pitch to neighboring pitches (for example, the note Sol can be adjusted to Fa or La as shown in Fig. 4) with a probability of HMCR × PAR (0 PAR 1) while the original pitch obtained in HM consideration is just kept with a probability of HMCR × (1-PAR).
where i x′ obtained in HM consideration and ) (k ) is a neighboring index used for discrete decision variables ( m has normally +1 or -1).

Violated harmony consideration:
Once the new harmony ) ,..., , ( is obtained using the above-mentioned three rules, it is then verified whether it violates problem constraints. Although the new harmony violates the constraints, it can be still kept in HM by taxing penalty, just as rule-violated harmony (for example, parallel fifth violation in Fig. 5) was still used in composition by the famous composers such as Ludwig van Beethoven. This technique is used for considering constraints in this study.
is better than the worst harmony in the HM, judged by objective function value, the new harmony is included in the HM and the existing worst harmony is excluded from the HM.

SCHOOL BUS ROUTING PROBLEM
The music-inspired HS algorithm is applied to school bus routing problem which is a multi-objective problem to minimize both the number of operating buses and the travel time of all buses, with two major constraints (bus capacity and time window). The study network to be optimized consists of one bus depot, one school, and ten bus stops as shown in Fig. 6. Each bus stop is demanded by certain number of commuting students, and travel time (in minutes) between two stops is specified in the figure.
While each decision variable i x stands for each bus stop (demand node) i , and has a value of the specific bus number k ( VS k ∈ ) that serves the bus stop i , the objective function of the problem is to minimize both the total number of operating buses and the moving time of the buses as described in the first and the second terms of Equation 6. The third and fourth terms represent penalty costs for the violation of bus capacity and time window, respectively. The fixed cost per each operating bus, fc is assumed as $100,000/bus; routing cost per moving time, rc is $105/min; shortest path ij sp (in minutes) between node i and node j is calculated by Floyd and Warshall's algorithm; connection status between node i and node j for bus k , k ij lk has 1 when bus k passes along from node i ); penalty cost for every bus capacity violation, 1 pc is $100,000 when any bus carries more than 45 students; and penalty cost for every time window violation, 2 pc is $100,000 when any bus operates more than 32 minutes to convey students to the school. Equation 7 represents the maximum number constraint for the operating buses, where the number of operating buses can be equal to or less than the number of candidate buses (= 4 buses in this study). Equation 8 shows that the number of boarding students in a bus is less than or equal to the bus capacity (= 45 students for each bus). Equation 9 shows that the travel time of a bus is less than or equal to the time window (= 32 minutes for each bus). Travel (in-vehicle) time in this study consists of both bus moving time and student boarding time, where each student's boarding time bt is 6 seconds.

COMPUTATION AND RESULTS
In order to apply the HS algorithm to the school bus routing problem, problem and algorithm parameters are specified: number of musical instruments (number of demand nodes) = 10; pitch range of each instrument (value range of each decision variable) = {bus 1, bus 2, bus 3, bus 4}; HMS = 10 ~ 100; HMCR = 0.3 ~ 0.95; PAR = 0.01 ~ 0.05; and stopping criterion = 1,000 improvisations. Next, harmonies (solution vectors) are randomly generated from the possible range as many as HMS. After that, a new harmony is improvised based on three rules (random selection, HM consideration, and pitch adjustment). The new harmony vector x′ is then put to the objective function to obtain total cost vtm , becomes 1, and the penalty cost for the time window violation is added to the objective function. If the total cost for the new harmony is better than the cost of any harmony stored in the HM, the new harmony is included in the HM, and the existing harmony with the worst cost is excluded from the HM. The new harmony is continuously improvised until the stopping criterion is satisfied.
When applied to the school bus routing problem, HS algorithm could find the global optimum ($307,980) only after 1,000 improvisations (or function evaluations), which was examined by the total enumeration for this 4 10 ( 1.05 × 10 6 ) combinatorial problem. It took 6.6 seconds on Intel 233MHz CPU to search 0.1% of total solution space (= 10 3 / 1.05 · 10 6 ). Table 1 shows the optimal route, number of commuting students, and travel time for the global optimal solution and near-optimal solution. "Do Nothing" in 5 th row of the table means that the bus 4 does not have to be operated because only three buses can serve all the students in the problem.
HS results were also examined by comparing with those of popular evolutionary technique, GA. In order to fairly compare HS with GA, the number of objective function evaluations and the number of computational runs are same in both algorithms: in HS, the number of function evaluations is 1,000 and the number of runs is 20 with different HMS, HMCR, and PAR; and in GA, the number of function evaluations is 1,000 (= population size × number of generations) and the number of runs is 20 with different population size, crossover rate, and mutation rate, recommended by Koumousis and Georgiou [21] .
After being applied to the test routing problem, both algorithms could find the global optimum solution. However, HS reached it twice while GA did it once out of 20 different runs. The average costs are $399,870 in HS, and $409,597 in GA, respectively. In addition, HS performed each run slightly faster than GA on the same machine.

CONCLUSION
A recently-developed music phenomenon-inspired algorithm, HS was introduced and modeled for solving the school bus routing problem. The objective of proposed HS model for the school bus routing is to minimize the total cost of multi-objective function which consists of bus operating cost, bus travel time, and penalties related with bus capacity and time window violations. HS model could find global optimum within far less function evaluations comparing with total enumeration. HS model also found better solution than GA in terms of number of reaching global optimum, average cost out of multiple runs, and computing time.
From these results, the HS algorithm appears to have a potential to be successfully applied to combinatorial problems in traffic engineering field. Especially, the presented bus routing model is expected to be applied to large-scale routing networks for the future study.