173  Advances in Production Engineering & Management ISSN 1854 ‐6250 Volu me 16 | Number 2 | June 2021 | pp 173–1 84 Journal ho me: a p em‐journal.or g https://doi.org /10.14743/apem2021.2.392 Original s cientif i c paper     Improved Genetic Algorithm (VNS‐GA) using polar  coordinate classification for workload balanced multiple  Traveling Salesman Problem (mTSP)  Wang, Y.D. a , Lu, X.C. b,* , Shen, J.R. c   a Beijing Jiaotong University, Shangyuan Village, Haidian District, Beijing, P.R. China   b Beijing Jiaotong University, Shangyuan Village, Haidian District, Beijing, P.R. China  c Beijing Capital Agribusiness & Food Group Co., Ltd., Xicheng District, Beijing, P.R. China       A B S T R A C T   A R T I C L E   I N F O The multiple tr a veling salesman problem (mTSP) is an extension of the travel‐ ing s a les m an p roblem ( TSP), which has wider applications i n rea l life than th e trav eli n g salesman p ro blem s uch as t ransportation a n d delivery, t ask alloc a ‐ tion, e t c. I n thi s paper, a n improved g enetic a lgorithm ( VNS‐GA ) that uses p o lar coordinat e c las s i f i cation to generate the i nit i al s ol utio n s is p r op os ed. It integrates the v ariable neighbourhood a lgori t h m t o solve the m u ltiple o bjec‐ tive o ptim izati o n of the m TSP with worklo ad b alance. Aiming t o workload b a l a n c e , t h e f i r s t d e s i g n o f t h i s p a p e r i s a b o u t g e n e r a t i n g i n i tial s olut ion s bas e d on the p olar c oordinate c la ssification. Then a d istance c omp a rison ins e rtion op erator i s des i gned a s a neighbourhood a ction f o r al locating p aths in a t argeted manner. F inall y , t h e neighbourho o d descent proces s in the v ari‐ a b l e n e i g h b o u r h o o d a l g o r i t h m i s f u s e d i n t o t h e g e n e t i c a l g o r i t h m for t h e e x p a n s i o n o f s e a r c h s p a c e . T h e i mproved algorithm i s tested o n the TSPLIB standard d ata set and compa r ed w ith other genetic algorit h ms. T he results show that the i mproved gene tic algorit h m c a n increase c omputati o n al e ffi‐ c i e n c y a n d o b t a i n a b e t t e r s o l u t i o n f o r w o r k l o a d b a l a n c e a n d t h is a lgorith m ha s w ild a p p lica t ions in r e a l l if e s u c h a s m u ltip le r ob ots ta s k allocatio n, scho o l bus routing problem and othe r optimization problems.   Keywords: Multiple traveling salesman pro blem (mTSP); Workload balance; Variable neighb o urhood search algorithm (VNS ) ; Genetic algori th m (GA); Polar coordinates; Classi fi cati o n *Corresponding author: xclu@bjtu . edu.c n (Lu, X.C.) Article history: Received 4 June 2021 Revised 15 June 2021 Accepted 17 June 2021 C o n t e n t f r o m t h i s w o r k m a y b e u s e d u n d e r t h e t e r m s o f t he Cr eat i v e C om mon s A t t r ib ut ion 4.0 In t er n at ion a l L i c e n c e ( C C B Y 4 . 0 ) . A n y f u r t h e r d i s t r i b u t i o n o f t h i s w o r k m u s t m a i n t a i n a t t r i b u t i o n t o t h e a u t h o r ( s ) a n d t h e t i t l e o f the work, journa l c i t a tion and DOI.      1. Introduction   Multiple T raveling S ales men Probl e m ( m T S P ) i s a n e x t e n s i o n o f t he clas s ic T raveling S ales man Proble m (T SP)[1 ] . T h e m T SP with workload b alance is different f ro m t h e basic TS P. I n som e scenarios, the b al anced distribution of w orklo a d is a n impor tan t consideration. F or e xample, school bus r outing probl em, rout e ar rangem ent p r oblem, p ost m an distribu tion problem, express delivery a rea division, e tc., can all be a bstracted as mTSP wit h balanc ed wo r kload. Omar et al . [2] proposed that mTSP c an b e used to s o lve vehicles, r obots, a nd U AVs routing problems. In s ingle‐ objective opt i mizati on o f m TSP, th e to tal c ost is the m ost co mm only used o ptimization objective. It c aus e s th e problem o f a n u n b a lanced d istribution of d is tance a m o n g t r a v e l i n g s a l e s m e n , w h i c h i s d i r e c t l y r e f l e c t e d i n t h e p a t h a s t h e d i s t a n c e t r a v e l e d b y o n e t r a v e l e r i s l o n g e r t h a n o t h ‐ e r t r a v e l e r s . T h i s s i t u a t i o n i s h o p e d t o b e i m p r o v e d i n s o m e i n dustries in reality bec a use the Wang, Lu, Shen    174  Advances in Production Engineering & Management 16(2) 2021 unbal a nced w orkload distribution res u lts in oper a tional difficu lties. Moreover, avoiding crossing paths with o ne a nother i s the most s i g nific a nt problem in prac t ical w ork. F or instance, i n school bus routing problem, if a certain r o u t e i s t o o l o n g d u r i n g p l a n n i n g s c h o o l b u s r o u t e , i t m a y c a u se the first arrival point to a rrive t oo e arly, which will affect students’ study and other perfor ‐ manc e. I n p o stmen ro utin g problem, a voiding cros s ing paths can improve the efficiency o f work. T h e r e f o r e , o u r m a i n f o c u s i s o n t h e m T S P w i t h w o r k l o a d b a l a n c e . T o solve this problem, we build a multi p le o bjective o ptimization mod e l and propose an i m proved g e n etic a l g orit hm th a t aims to mi ni mize the tota l cost and bal a ncin g the w o rkload in t his paper. 2. Literature review  R e s e a r c h o n T S P a n d m T S P i s o f l o n g h i s t o r i c a l m a k i n g . T h e o r i g i n a l s o l u t i o n t o t h i s t y p e o f problem mainly u sed acc u rate a lgorithms. A li et al. [ 3 ] u s e d t h e b r a n c h a n d b o u n d m e t h o d t o solve t h e symmetrical m TSP wit h 5 9 cities. Gavis h et al . [4] d e scribed that l imiting the low e r bound o f th e b ranch i m proved the b ranch‐and‐bound meth od, which c a n s o l v e t h e m T S P w i t h 1 0 t r a v e l i n g s a l e s m e n a n d 1 0 0 c i t i e s . H o w e v e r , a s t h e s c a l e o f t h e p r o b l e m e x p a n d s , t h e s c a l e o f the solution s pace a lso increas es e xponenti a lly, makin g i t diff icult to u se a ccurate a l g orithms. There f or e, h euristic algo r ithms have b ecom e th e choice o f mo re and mor e r esearchers. Russell [5] was the first one to apply the heuristic algorithm t o s o l v e t h e m T S P , w h i c h t r a n s ‐ f o r m e d m T S P i n t o T S P a n d t h e n s o l v e d i t w i t h t h e L i n ‐ K e r n i g h a n algorith m. I t also o b tained a more a ccurate a pproximate solution . For solving the mTSP w ith t he s hort est total dist ance a nd the shortest s ub‐path, C arter et al . [ 6 ] designed a two‐stage c hromoso m e coding s cheme and d e s i g n e d a g e n e t i c o p e r a t o r t o i m p r o v e t h e s o l u t i o n s p e e d . S i n g h et al. [ 7 ] u s e d t h e i n v a d i n g w e e d a l g o r i t h m t o o b t a i n a c o m p a r a t i v e o p t i m a l s o l u t i o n . Z h o u a nd L i [ 8 ] used a m utat ion oper‐ a t o r c o m b i n e d w i t h a 2 ‐ o p t s e a r c h o p e r a t o r t o i m p r o v e t h e g e n e t ic a lgorithm, avoi ding th e premat ure p h eno m en on of g en etic a l g orithm, and used s imul ation method s to v erify t he effec‐ tiven e ss of t he a lgorith m . Xu et al . [9] also c o m bin e d the 2‐opt operator w it h a genetic alg o rithm t o s o l v e t h e m T S P . H u et al. [ 1 0 ] p r o p o s e d a n i m p r o v e d g e n e tic algorith m that i ntegrates t he reproduction mechanis m of t he w eed algorithm and the locally o p timized mutation o perator to s o l v e t h e m T S P w i t h b a l a n c e d w o r k l o a d . G u o [ 1 1 ] u s e d a g e n e t i c algorith m to a nalyze a nd v eri f y the design o f chromosome c oding scheme f or t he g enetic a lgor ithm t o s o l v e m T S P . L u et al. [ 12] adopted a tw o‐stage al gor i thm to solve m TSP while avoidin g th e path‐crossing probl em between traveling sal e smen. Bostanci et al. [ 1 3 ] u s e d t h e i n t e g e r l i n e a r p r o g r a m m i n g m e t h o d t o s o l v e t h e problem o f f inding t he s hortest ro ut e o n n etw o rks with G IS. Thi s study provided a n ew m ethod t o s o l v e t h e o p t i m i z a t i o n p r o b l e m s , h o w e v e r , o n l y T S P w a s s o l v e d b y th is meth od, problems more c ompli c ated such a s mTSP ar e d i fficult to fi n d the opti m a l solution. R e c e n t l y , w i t h t h e r i s e o f v a r i o u s h e u r i s t i c a l g o r i t h m s a n d m a c hine l earn ing, i n addition to genetic algor ithms, ant colony o ptimization al gorithm (ACO) [14], a rtificial bee colony a lgorithm (ABC)[15], t abu search a lgorithm (T S )[1 6 ] a nd o the r he u ristic a l g orithms have a lso been tried t o s o l v e m T S P . S o n g [ 1 7 ] u s e d t h e s i m u l a t e d a n n e a l i n g a l g o r i t h m ( S A) to solve the problem of three traveli ng sales m en i n 40 0 cities, but the algorithm took a long time . Hu [1 8] constructe d a n a r c h i t e c t u r e c o m p o s e d o f a s h a r e d g r a p h n e u r a l n e t w o r k a n d a d i stributed strategy n etwork to gen e rat e a n approxi m at e optimal s o lution for mTSP a nd u sed rein forcement learning t o train the model, a nd this met h od s hows g ood results in large‐scale e x amples. Justus [ 19] abstracted the problems i n th e p a th o f staff guiding t o urists d uring th e M ecca p ilgri m age as m ulti ple obj e c‐ t i v e m T S P w i t h t i m e w i n d o w a n d p r o p o s e d a n i n t e r a c t i v e m e t h o d f or p roviding a s oluti on. A im‐ ing the mTS P u nder d ifferent g oals, Liu et al . [20] p roposed the For e stTraversal al gorithm and the Retr ace algorith m t o s olve the m aximum‐minimum m TSP and the m T S P t h a t m i n i m i z e s travel costs, respectively. Genetic al gorithm (G A) i s a h e uristic al gorith m t h at is often u sed wh en s olving m TSP . T hi s method m ai nly focuses on codin g d esign [5] an d how t o a void g en etic a lgorithms prematurely falling into l o cal optimal ity. Th e variable n eigh b o urhood s ear c h algorithm (VNS) has t h e ability to e xp and th e se arch r an ge b y s y ste m atically c h a ngi n g th e n e igh bourho o d s tructure. It g u a ran‐ t e e s t h e l o c a l s e a r c h c a p a b i l i t y w h i l e e n s u r i n g s o l u t i o n s f o r h aving good d iversity [ 21]. T here ‐ Improved Genetic Algorithm (VNS‐GA) using polar coordinate classification for workload balanced multiple Traveling …   Advances in Production Engineering & Management 16(2) 2021  175 fore, this p aper d esigns a g en etic a lgorithm b ase d on p o lar c o o rdinates t o quickly generate h i g h‐ quality i n itial solutions, and com b i n e s V N S t o r e t a i n t h e c h a r a c t eristics o f neighbourhood d iver‐ sity, designs different n ei ghbo urhood a ctions to i m prove the al gorithm's search a bility to solve mTSP with workload balance. 3. Mathematical model and problem solving  This p aper p roposes an i mproved genetic al gorithm, w hich i s imp roved by f using the variable neigh b ourh o o d al gorith m ( V ari a ble nei g hb ourh ood s earch g en etic algorithm), and VNS‐GA is referred t oo, as an abbr eviation in t h e text. 3.1 Problem description and mathematical model  T h e m T S P s t u d i e d i n t h i s p a p e r c a n b e d e s c r i b e d a s g i v e n n c i t i e s , m t r a v e l i n g s a l e s m e n a n d t h e distance m atrix betw een nodes D ij = (d ij ) n·n , t r avel er s s t ar t fr o m the same node (sourc e node), visit a certai n nu mber o f nodes a nd t h e n retur n t o the source no d e . A l l n o d e s a r e r e q u i r e d t o b e a c c e s s e d , a n d t h e r e s t o f t h e n o d e s e x c e p t t h e s o u r c e n o d e c a n only b e a ccessed once. T he a im is to f ind a H a miltonian c ir cuit that makes total distance s hortes t and worklo ad b al ance. Assuming the problem is s ymmetric al m TSP, t hat is, d ij = d ji , t he multiple objective op timizati on m odel can be described as follows: 1 , 2 ,…, (1) m a x m i n , , ∈ , ∈ ( 2) ∀ 0 , 1 ,…, ; 1 , 2 ,…, (3) ∀ 0 , 1 ,…, ; 1 , 2 ,…, (4) 1, 1 , 2, … , ,0 ( 5 ) ∈ (6) | ∉ ∈ 1 ,⊂ 0 , 1,2, … , , ∀ 1,2, … , (7) 1 traveler wen t through a rc 0 otherwise (8) 1 traveler v isited node 0 otherwise (9) Where V = {0,1,2,...,n} is t he node set, n i s t h e n u m b e r o f n o d e s ; T = {1,2,...,m} i s t h e t r a v e l e r s set, m is the traveling sal e sman n umber; A = {( i, j, k)|i, j = 0,1,2 , …,m, i ≠j , ∀k = 1,2,…,m } is t he a rc s e t , w h i c h r ep r e s e n t s t he s e t of r ou t e s t h a t th e t r a ve l e r s ma y pass; D = {d ij |i, j ∈ A } i s t h e d i s t a n c e matrix, r e presentin g th e d istance fr o m n ode i to n o de j . E q s . 1 a n d 2 a r e t h e o b j e c t i v e s o f t h e s e f u n c t i o n s . E q . 1 m i n i m i z es t h e t otal p ath distanc e ; Eq. 2 minimi zes the differen ce b etwee n t he l o n gest p ath and t h e shor t e s t p a t h i n m p a t h s . E q . 3 means that f or a ny d own s tream n o de, only on e up stream n ode is a llowed for its connec tion, a nd Eq. 4 means that in each u pstream no de, only one d ownstream nod e is a llowed to b e c o nnected; Eq. 5 me ans that e xcept for the source node, a ll other nodes ar e o n ly v isite d once, a nd a ll travel‐ ers start fro m t he s ourc e node. I n E q . 6, S i s t h e b r a n c h e l i m i n a t i o n c o n s t r a i n t , a n d E q . 7 i s a n expression o f S [ 2 2 ] , w h i c h m e a n s t h a t f o r a n y n o d e i n t h e t r a v e l e r s ’ p a t h s , a ny o f its subsets must b e co n n ected to the other su b set s in the solut i on. Wang, Lu, Shen    176  Advances in Production Engineering & Management 16(2) 2021 3.2 Fitness function  In m ultiple objective op timizati on p r oblems, the goals o f ten co nflict w ith each o ther. The most common method is the g en eratio n me thod ( inc l uding weighting met ho d, c onstraint method, e t c . ) , i n t e r a c t i v e m e t h o d , a n d h y b r i d m e t h o d t o s o l v e m u l t i p l e o b je ct i ve s. T h e we ig h t in g me t h o d assigns different wei g hts to m ultiple goals which can convert t he multiple objective opt i mizati on problem in to a s ingle‐ obj e ctive opti mization pr obl e m and mak e s the simpl e a nd e asy i m plem en‐ tation of the p roblem. Therefore, th e w ei ghtin g m ethod is u sed in t his paper to t ransform t he two obj e ctives optimi z ati on probl e m into a sin gle‐ objective opt imizati on pr oblem. (10) Where ∑∑ , 1 , 2 ,…, . w 1 , w 2 are the w e ight coefficients, which b elong to hyperpar am eters and th e valu es o f c o efficients v ary with t he d at a s e t . I t c a n b e d e t e r m i n e d b y expert evaluation method o r comparison method, etc. Therefor e, t he fit ness fu nction in the gen etic al g or ithm is: 1/ (11) I t c a n b e seen f rom Eq. 11 t ha t the g re a te r the f i tne ss, t he str o n g e r t h e i n d i v i d u a l a d a p t a b i l i t y , and th e b e tt er the repres e nted f e a sibl e solution. 3.3 Chromosome coding  I n t h e g e n e t i c a l g o r i t h m , t h e d e s i g n o f t h e c h r o m o s o m e e n c o d i n g n eeds to refl ect th e genetic characteristi c s of th e ind ividual and s h ould b e easy to operate . To e xpress the pat h m o r e e ffec‐ tively, a o n e ‐ piece decim a l codin g s c h em e is a do pted i n this p a per. A c hromoso m e sh ould in‐ clude a t least one node except th e s o u rce node, a s shown bel o w: T a k e t e n n o d e s a n d t h r e e t r a v ‐ elers as a n ex ample, u se 0 , 1, 2 , 3, 4 , 5, 6 , 7, 8 , 9 as it res pecti v ely repr esen ts ten nod es, of w hich 0 i s t h e s t a r t i n g n o d e o f t h e t h r e e t r a v e l e r s . T h e n t h e p a t h o f a t raveler can be e xpressed as: 0‐1‐ 2‐3‐4‐ 5‐6‐7‐ 8‐9. To r epresent multiple travelers, virtu al nod es 1 0 a nd 1 1 a re i n serted to repr esent the starting node o f t h e r oute, which is a lso th e s o urce n ode. T hen, the n ew c hromoso m e c a n be f o r med as: 0‐1‐2‐ 3‐4‐5‐ 10‐6‐ 7‐1 1 ‐8‐ 9 . T a k i n g v i r t u a l n o d e s 1 0 a n d 1 1 a s the path d ividing nodes, then t h e t h r e e p a t h s c a n b e e x ‐ pressed as: ( 1 ) 0‐ 1‐2‐ 3‐4 ‐ 5‐0, (2) 0‐6‐ 7 ‐0, (3) 0‐ 8‐9‐0. 3.4 Polar coordinate classification method initialize population  The qu ality of t he i niti al p opulation p l ays an i mpo r tant r ole i n the intelligent g roup o ptimization a l g o r i t h m . F o r m T S P w i t h a s i n g l e s t a r t i n g n o d e a n d c l o s e d c y c l e, i t is a c ommon opti mizati on method t o fi rstly use the clustering method t o optimize g rouping , a n d t h e n u s e a h e u r i s t i c a l g o ‐ r i t h m t o o p t i m i z e t h e o r d e r w i t h i n t h e g r o u p . I n 2 0 1 9 , M i n et al. [ 2 3 ] p r o p o s e d a m e t h o d u s i n g MRISA ( m ult i ‐restart‐iteration sweepi ng) and tabu s earch to solv e v e h i c l e r o u t i n g p r o b l e m . T h i s two‐stage al gorithm also c an b e d e scribed as clustering b y d ist ance in polar coordinates and o p t i m i z a t i o n b y t a b u s e a r c h . T h e n i n t h e s a m e y e a r , t h e y i m p r o v e d t h i s a l g o r i t h m b y a d d i n g a n a d j u s t o p e r a t i o n n a m e d p u t i n & p u t o u t . T h i s i m p r o v e d a l g o r i t h m has thr ee stages: cl uster the customer p o i nts, u se p ut i n & put out operatio ns t o adjust the load d em and and use t a b u s earch t o s o l v e t h e T S P [ 2 4 ] . T r a d i t i o n a l c l u s t e r i n g a l g o r i t h m s s u c h a s K‐means cl ustering o r fuzzy C ‐ means clustering and other clustering methods based on centroid s (although they can q uickly cluster the s et o f nod es that ne ed to b e v isited) can’t well r e present the starting poi nt a s th e origin wher e e ach traveler’s s ub‐p ath is d istributed a round the s ame startin g p oint. Th erefore, a classification method bas e d on po lar c oordinates is proposed in this paper. The specific o peratio n p r o cess of the p olar c oordi nate classifi c a tion algorithm to g enerate t h e i n i t i a l p o p u l a t i o n i s a s f o l l o w s . F i r s t , e s t a b l i s h a p o l a r c o o r dinate s ystem with the s tar ting node a s t h e o r i g i n a n d m a p n o d e ‐ s e t (n n o d es) into the p olar c oordinate s ystem then s ort each node by t he a ngle i n th e pol a r coordi nate s ystem. T hen calculate t he two nodes w ith the furt hest a n‐ gular distances and select one o f them to connec t with the o r i g in t o est a b l ish a new p o lar coor‐ Improved Genetic Algorithm (VNS‐GA) using polar coordinate classification for workload balanced multiple Traveling …   Advances in Production Engineering & Management 16(2) 2021  177 d i n a t e s y s t e m a x i s . A f t e r t h a t u s e t h e n e w a x i s a s t h e s t a r t i n g p osition to r otate clock wise, and use the trav eler n u m b e r (m t r a v e l e r s ) t o e q u a l l y d i v i d e t h e n u m b e r o f n o d e s t h a t t h e t r a v ele r s n e e d t o v i s i t . T h e n o d e s i n t h e a r e a s w e p t d u r i n g t h e r o t a t i o n are the n o de ‐set t o be v is ited b y a traveler, and m i n i t i a l p a t h s ( t h a t i s , e v e r y t r a v e l e r r o u g h l y n e e d s t o v i s i t n/m nodes) c onstitu te an i nitial f easible solution. Finally, us ing chr o mos o me f r a gmen t flippin g a s nei g hborh ood s hak‐ ing to g en er ate e nou gh c hromoso m es a s the initial population. T he p rocess o f the pol a r coordi‐ nate classific a tion algorithm is as shown in Al g orithm1: Algorithm1 ：Polar C oor dinate i nitial ize pop ulation I n p u t : T h e n u m b e r o f p o p u l a t i o n p , t h e nu mb er of sales m en n , cartesian coordinates of th e cities (x i ,y i ) and the num b er of cities m Output: Initialized population k = 1 ( θ i , d i ): chan g e c a rtesian coordinates (x i , y i ) to pol a r coordinat e s and s o rt t he cities by t he an gl e θ i max ( θ, d ): c alculate the maxim um angle differenc e am on g t h e cities and build new polar c oordinate by joining the origin city a n d one of ma xim u m a n g le difference cities s 0 : gener a te a solutio n b y spinnin g clockwise and each sales m an v isits n/m cities repeat s : neighborhood shaki ng s 0 k = k + 1 until k > p end 3.5 Neighbourhood actions in VNS  Cross‐recombinatio n and mut a tion o perators a re o ften u sed in g e netic al gorithm to g enerate new p o pulat i ons. To e nable th e algor i thm in f i n ding a m or e effi cient com b ination th at m ak es t h e travelers ’ w orkload ro u g hly b a lanced, the distan ce c o m paris on i nsertion o per a tor co operates with cross‐recombination (neighbo u r hood a cti o n 2) a nd m utation op er ater ( nei g hb ourhood action 3) is d e signed in th is pap er to p r oduce the n e xt g e ner at i on more effi c iently. Neighb ourh ood a cti on 1: D istance co mparison o p e rator to f in d a bett er s o l ution for making the workload b alance f aster. T he n ei ghbo urhood a ction adopted in t h i s p a p e r i s t o r a n d o m l y s e l e c t a n o d e f r om t h e s ub ‐ p a t h w i t h th e l on g e s t d is t a n c e a n d i nsert it i nt o the sub‐path w ith the shortest d istance. A s shown in Fi g . 1 : A s s u m e t h a t p a t h 0 ‐ 1 ‐ 2 ‐ 3 ‐4‐5‐ 0 h as the l ongest d istance amo n g all su b‐paths, a nd 0 ‐ 6 ‐7‐0 is t h e p a th w ith the shortest d i s t a n c e a m o n g a l l s u b ‐ p a t h s ; I n t h e d o m a i n a c t i o n , r a n d o m l y s e l e c t a n o d e i n t h e l o n g e s t p a t h ( s e l e c t n o d e 3 ) a n d i n s e r t i t i n t o the sh ortest path to f orm a new chromosome. If t he f itness of t he n ew c h r omoso m e i n creases after inserti n g t h e s e lec t ed node, t h e n ew c hro m osom e wil l t hen b e retained. Otherwise, t he paternal c hromosome will be r etai ned and this inferior solution w i l l b e r e c o r d e d t o p r e v e n t d u ‐ plication. If the f itn e ss of the n ew c hromoso m e increases a fter i nserting the selected node, s uggests t h a t t h e n e w o f f s p r i n g c h r o m o s o m e s generated after the neighbou rhood action are: 0 ‐ 1 ‐2‐4‐ 5 ‐ 10‐6‐ 3‐7‐ 11‐ 8‐9. Fig. 1 Schematic of gene insertion in neighbour h ood actio n 1 Wang, Lu, Shen    178  Advances in Production Engineering & Management 16(2) 2021 To a void r epeated calculations, only the differ e nce between the d i s t a n c e o f t h e s e q u e n c e i s calculated each time in the c a lculation process. A ssuming t hat the selected node is t he node at the i ‐ th p osit ion in the n ‐ t h p ath, a nd t he i nserti on position is t he n ode at t he j ‐ th p ositi on in the m ‐t h p a t h , t h e n a c t u a l d i s t a n c e t o b e c a l c u l a t e d i s d = d n – (d i ‐ 1,i + d i,i+1 ) + d i ‐ 1,i +1 – (d m + d j ‐ 1,j + d j,j+1 – d j ‐ 1,j +1 ) . Where d n , d m a r e t h e p a t h l e n g t h b e f o r e p a t h n a n d m t o perform th e n e igh b our h ood a c‐ tion, i – 1 and i + 1 are th e previo us a nd next nod es of i i n p a t h n , the same a s j – 1, and j + 1 is th e previous a n d n ext nodes o f j i n p a t h m , respectively. d ≥ 0 m e a n s t h a t t h e fitness is not increased after performing the nei ghbo urh ood a ction, oth er wise, it also m eans th a t t he fit n ess is i ncreased, and th e curr ent soluti on is better than the parent solution. Neighb ourh ood a ction 2: l ocal crossover oper ator . Neighbourh ood a c t i o n 1 c a n m a k e t h e o f f ‐ s p r i n g d e v e l o p i n t h e d i r e c t i o n o f w o r k l o a d b a l a n c e , b u t i t i s easy t o fall i nto the loc a l optimum. Therefor e, a n improved crossover op erator to avoid the offsp r in g from f al ling i nto the local op‐ timum is proposed in this paper. T o a v o i d r e p e a t e d f r a g m e n t s o r m i s s i n g f r a g m e n t s i n t h e o f f s p r i ng t hat ar e generated afte r the crossover and for quickly completing t he crossover, we firs t randomly s elect a sample g ene f r a g m e n t f r o m t h e p a r e n t 2 a n d i n s e r t s i t i n t o t h e o f f s p r i n g a t th e o riginal posit i on. T hen t r a v e r s e s pa r e n t 1 t o f ind t h e g en e s th a t a r e d if f e re n t f r o m t h e sa m ple g ene f ra g m e nt a nd i nse rt‐ e d t h e m i n t o t h e o f f s p r i n g i n s e q u e n c e , t h e g e n e s t h a t a r e t h e same a s the sample g ene fragme nt is skipped. The crossov e r process is shown in Fi g. 2 . Neighb ourh ood a ctio n 3: m u tatio n o p erator. Jud g ing fro m t h e c ha r a cteristics o f th e tr aveling salesman p roblem, ther e should be n o c r o s s e d g e s i n t h e r o u t e o f th e opti mal solution. B ased on t h i s f e a t u r e , a s s h o w n i n F i g . 3 (a i s b e fore m ut ati on, b i s aft e r mutati on), the destruction & re‐ construction m ethod is u sed to m utate individuals to reduce cro ssover between p aths. Suppose the p a th s et i s S = {s 1 , s 2 ,...,s i ,. ..,s k ,. ..,s n }, t hen sort t he d istance ma trix by dista nce t o ge t so rte d n o d e pairs and r a ndomly s electing a node p i n a p a t h s i a nd i nterch ange i t with t he n ear e st n eigh bour node m i n t h e nearest nei g hbo u r rout e s k t o r e d u c e t h e i n t e r s e c t i o n b e t w e e n p a t h s . I f t h e f i t n e s s i n c r e a s e s , u p d a t e t h e c u r r e n t s o l u t i o n ; o t h e r w i s e , i t w i l l n o t u p date a nd r ecord the nod e pair to node p a ir distance matrix for a voiding repeated exchanges. Fig. 2 Cro ss pro cess di agr am Fig. 3 Mutation operator diagram 3.6 Generate offspring   We i ntroduc e d the idea o f searching the neighbourhood s tructure a l t e r n a t e l y i n V N S t o e x p a n d the search s pace o f GA, along w ith t h e elite r e te ntion strat e gy w h i c h i s u s e d f o r r e t a i n i n g t h e dominant in d ividuals to i m prove the conver genc e speed. The specific process of g enerating the next popul ation is f irstl y f o r c a l c u l a t i n g t h e f i t n e s s o f each c hromo s ome in th e c u r r e n t p o p u l a t i o n a n d r e c o r d t h e o p t i m a l i n d i v i d u a l . T h e n , us e the r ou‐ Improved Genetic Algorithm (VNS‐GA) using polar coordinate classification for workload balanced multiple Traveling …   Advances in Production Engineering & Management 16(2) 2021  179 lette to selec t s epar ate c h romoso mes and search in t he nei g hb ourhood structure N k (k = 1 , 2 ,...,n ) t o f i n d a b e t t e r s o l u t i o n . T h e b e t t e r s o l u t i o n w i l l b e a d d e d t o the p rogeny p opulation and k w i l l b e s e t t o 1 . I f a n y b e t t e r s o l u t i o n c a n n o t b e f o u n d , s k i p t o t h e next n ei gh bourhood, and repeat the process of sel ection a nd n ei gh bo urhood desc e nt u ntil th e si ze of th e p o pulation reaches max‐ i m u m . A t t h i s t i m e , i n d i v i d u a l s w i t h g r e a t e r f i t n e ss c a n b e se l ecte d more i n the se le ction proce ss a n d t h e i n d i v i d u a l ’ s n e i g h b o u r h o o d s p a c e c a n a l s o b e f u l l y s e a r ched. Th e pseudo‐code of th e generating next p o pulati on algo rith m is shown i n Algorithm 2 : Algorithm2: Variable neighbourhood descent to genera te n e x t generation Input: origi n population S 0 Output: new population S s 0 : optimal solution of origin pop u lati on k =1 add s 0 to S S n = 1 repeat s : roulett e a solution fro m S 0 repeat s’ : fi nd best n e ighb our s’ if f(s’) > f(s) then add s’ to S k = 1 else k = k + 1 until k = k max until S n = S max end 3.7 Algorithm flow  The VNS‐GA algorithm p r o cess is as below: 1) U se p olar c oordinat e classification a l g orithm t o gen e rate i n itial pop u lation S 0 ; def i ne s n neigh b ourh o o ds and den oted as N k (k = 1 ,2 ,...,n ). 2 ) V e r i f y t h e f e a s i b i l i t y o f t h e i n i t i a l p o p u l a t i o n , a n d u s e n e i g h b ourhood s h a king t o eli m inat e infeasible sol utions. 3 ) C a l c u l a t e t h e f i t n e s s o f i n d i v i d u a l s i n t h e p o p u l a t i o n , a n d s e l e c t t h e c u r r e n t o p t i m a l s o l u ‐ tion s. 4 ) T h e v a r i a b l e n e i g h b o u r h o o d d e s c e n t p r o c e s s : s e a r c h o p t i m a l s olutio n s i n nei g hb o u rhood structure N 1 , if i t f i n ds be tt e r so lu t io n s' than s , then let s = s'. 5) I f no b et t e r solution c a n b e f ound i n the nei g h b ourhood s t r u cture N 1 , t hen let k + 1 and go to step 4. 6) Add the c u rrent opti m al solution to the n e x t‐gen e ration popu lation; 7) R epe a t st eps 4‐6 until the n u m b er o f i n dividuals in t he p op u lation r e ac hes the maxi mum value. 8 ) R e p e a t s t e p 7 u n t i l t h e m a x i m u m n u m b e r o f i t e r a t i o n s i s r e a c h e d , t h e n o u t p u t t h e c u r r e n t optimal sol u tion. 4. Results and discussion  To v erify th e perform a n c e of V NS‐G A, f irst, we q uoted th e data o f t h e C h i n e s e T r a v e l i n g S a l e s ‐ man Pr oblem (CTSP) i n literature [ 25]. U se P yth on code t o t e st it on a co mputer w ith the oper at‐ i n g s y s t e m W i n d o w s 1 0 , I 5 ‐ 8 2 5 0 , 8 G R A M . T h e a l g o r i t h m r u n s t e n t i m es i ndepende ntly o n each c a s e s e t . F o r t h e s a k e o f f a i r n e s s , t h e s a m e g e n e t i c a l g o r i t h m parameters a s in literature [25] a re Wang, Lu, Shen    180  Advances in Production Engineering & Management 16(2) 2021 u s e d , a n d t h e p o p u l a t i o n s i z e i s s e t t o 5 0 , t h e h e r e d i t y i s 1 0 0 0 generations, the crossover proba‐ b i l i t y i s 0 . 8 , a n d t h e m u t a t i o n probability is 0 .15. Because th e results obtained u sing SA in litera‐ ture [25] are wrong. According to th e d etailed p a t h dia gram o bt ain e d b y t he si m ulat ed a nne alin g algorithm in the literatur e, the r esults are revised a s shown i n Table 1. I n t h e c o m p a r i s o n o f t h e t h r e e a l g o r i t h m s , t h e t r a d i t i o n a l g e n e tic algorithm perfor ms the w o r s t a m o n g t h e t h r e e a l g o r i t h m s b e c a u s e o f t h e p r o b l e m e a s y f a lling into the local optimal so‐ lution. After u sing t he v ariable n e igh b ourhood a lgorithm t o i m p rove t he g en etic a lgo r ithm, the effect is pro m oted. C o m p ared w ith t he det a iled p ath diagram given in [2 5] ( F ig . 4 ) , i t c a n be s ta t ‐ ed th a t t he i m proved g enetic a lgorit hm is consis tent w ith t he n odes s et a ssigned b y th e sim u lat‐ e d a n n e a l i n g a l g o r i t h m i n t h e l i t e r a t u r e . B u t t h e i m p r o v e d g e n e tic algorithm has better p er for‐ manc e th an t he g en etic al g orithm on a single path, s o, the tot a l obtained path is shorter. In a ddition, t his paper sel e cts three d a ta s ets: e il51, kroA100, a n d k r o B 1 5 0 i n T S P L I B t o t e s t the per f ormance of t he a lgorithm . Th e performan c e of the a lgorithm u n der the conditi ons of 3, 5 , 10, a nd 2 0 t r avelers is t ested agains t the scal e o f t hree d at a sets, and th e p a rameter settings used in each d ata set are shown in Ta b le 2 throu gh e xperim en ts. We c omp a re d the m e asu r ed d ata se t results with the t est res u lts o f ot her algorith ms ( data f r o m l i t e r a t u r e [ 1 4 ] a n d [ 2 6 ] ) . T h e r e s u l t s a r e s h o w n i n T a b l e 3, w here n i s t h e n u m b e r o f n o d e s and m i s the numb er o f tr avel ers. G A1C is s ingle chromoso me c oding g enetic a lgorithm, GA2C is the two‐chr o mosome c oding geneti c algorithm, G A2PC is two‐segme nt c h r omoso m e c o ding g e‐ netic algorithm, GGA‐SS is s teady‐state grouping g enetic a lgorit h m , T C X i s i m p r o v e d t w o ‐ s e g m e n t c h r o m o s o m e c o d i n g g e n e t i c a l g o r i t h m , a n d R G A i s b e l t G e netic algorithm with r epro‐ d u c t i o n m e c h a n i s m , R L G A i s a g e n e t i c a l g o r i t h m t h a t i n t e g r a t e s the reproduction mech anism o f t h e w e e d a l g o r i t h m a n d l o c a l o p t i m i z a t i o n , I W O i s t h e i n v a s i v e weed a l g orithm, and VNS‐GA is the impr ove d gen etic al g orit hm prop o sed in this paper. Table 1 Comparison results between VNS‐G A and other algor ithms (CTSP da ta set) Algorithm T otal distance(km) Sub‐path distance(k m) V a riance GA 20225 6106 643261 6477 7643 SA Origin data Modified data Origin data Modified data Modified variance 17731 17404 5527 5527 1116314 5616 4909 6578 6968 VNS‐GA 1 7153 6968 1361290 4658 5527 Fig. 4 The specific solution p ath of VNS‐GA and SA Table 2 Parameter setting of different examples Standard data s e t Number o f travel‐ ers Population size Cross rate M utation rate eil51 3, 5, 10 60 0.8 0 .2 kroA100 3, 5, 10, 20 80 0.8 0 .2 kroB150 3 , 5, 10, 20 80 0.8 0 .2 Improved Genetic Algorithm (VNS‐GA) using polar coordinate classification for workload balanced multiple Traveling …   Advances in Production Engineering & Management 16(2) 2021  181 Table 3 Comparison of the longest path (VNS‐ G A with other g enetic algo rithms) Data n m GA1C G A2C GA2P C GGA‐S S TCX RGA RLGA I W O V N S ‐GA D (%) 51 51 3 234 275 203 161 203 174 167 160 165 3.13 51 5 173 220 164 119 154 125 118 118 121 2.54 51 10 140 165 123 112 113 122 112 112 112 0.00 100 100 3 14722 16229 13556 8542 12726 10231 10115 8509 8613 1.22 100 5 11193 11606 10589 6852 10086 7895 7812 6767 6445 4.76 100 10 9960 10200 9463 6370 7064 6234 6225 6358 5764 7.41 100 20 9235 9470 8388 6359 6402 6233 6211 6358 5395 13.14 150 150 3 19875 21067 19687 13268 18019 14886 14629 13168 10878 17.39 150 5 15229 15450 14748 8660 12619 8998 8927 8479 7711 9.06 150 10 12154 12382 11158 5875 8054 5723 5613 5594 5937 6.13 150 20 10206 10338 10044 5252 5673 5372 5251 5246 5750 9.61 The d a ta g iv en i n Table 3 are th e v a lu es o f th e lon g est sub‐p a t h. Under the t hree dat a sets and different t raveler numbers circumstances, the res u lts of V NS‐ GA a r e f a r s u p e r i o r t o t h e G A 1 C , GA2C, and G A 2PC al gorit hms. The co mparison wi th t he GGA‐S S algo rithm shows that o nly in the c a s e o f 3 t r a v e l e r s i n t h e k r o A 1 0 0 d a t a s e t , t h e r e s u l t s a r e s l ightly w orse a nd the r emai ning r e‐ sults are better than the results of G GA‐SS. W hen compared w ith t he T CX a lgorithm, the results a r e s l i g h t l y w o r s e o n l y i n t h e c a s e o f 2 0 t r a v e l e r s i n t h e k r o B 150 data s et, and the other results are bett er t h a n the TCX algorithm. I n the c a se o f t h e kr oB150 d ata s e t of 1 0 and 20 t ra velers, the results are sl ightly worse t han RGA and RLGA a lgorithms. We use D t o v i sua lly dem onstra t e the improvemen t or l ack of r esults w her e D = ( best_result_amon g_ other_geneti c_algorithms – V NS‐ GA_result)/best_result_amon g _oth er_gen etic_algo r ithms × 1 0 0 % . Thus, D > 0 m e a n s t h e r e s u l t i s i m p r o v e d a n d i s t h e b e s t r e s u lt a mong a ll genetic algorithms w h i l e D < 0 m e a n s t h e r e s u l t i s w o r s e t h a n t h e b e s t r e s u l t . T o v i s u a l i z e h o w t h e p e r f o r m a n c e o f t h e a l g o r i t h m c h a n g e s w i t h t h e numb er o f tr avel ers increase, w e intro d uce the par a me ter . is the av erage o f d iffere nt n u m ‐ bers o f trav elers (i is th e numb er o f tr avel ers and i = 3, 5, 10, 2 0). As c an b e see n f rom Fi g. 5, with the increas e of t ravelers’ number , the performance of the a lgori t h m d e c r e a s e , b u t i n g e n e r a l , t h e perform a nce of the algorithm is i m pro v ed co m pare d with other g e netic al gor i thms. Fig. 5 VNS‐GA performance variatio n wi th the n umber of trav elers Table 4 Comparison of total pa th (VNS‐G A with other genetic a lgorithms) Data n m GA1C G A2C GA2P C GGA‐S S TCX IW O VN S‐GA 51 51 3 529 570 543 449 492 448 495 51 5 564 627 586 479 519 478 552 51 10 801 879 723 584 670 583 1100 100 100 3 27036 30972 26653 22051 26130 21941 25524 100 5 29753 44062 30408 23678 28612 23319 31222 100 10 36890 65116 31227 28488 30988 27072 49606 100 20 62471 95568 54700 40892 44686 38357 101235 150 150 3 46111 48108 47418 38434 44674 38055 31960 150 5 49443 51101 49947 39962 47811 38881 35754 150 10 59341 64893 54958 44274 51326 42462 55505 150 20 94291 100037 73934 56412 62400 53612 106917 Wang, Lu, Shen    182  Advances in Production Engineering & Management 16(2) 2021 The data i n Table 4 are t h e v a lues o f the tot a l p a t h . The b a lan ce d egr e e ne eds to b e c o mpared with the lon gest su b ‐pat h valu e and t h e tot a l p a th valu e . Altho ugh on s o me data sets, th e value of the longest path i n th e s o lution obt a i n ed b y th e VNS‐GA a lgorit hm i s slight ly h igher th an o th er algorith ms, l iterature [ 14]is not m u l t i p l e o b j e c t i v e o p t i m i z a t i on and the objective of R GA a nd R L G A a l g o r i t h m s i n t h e l i t e r a t u r e [ 2 6 ] i s t o m a k e t h e l o n g e s t s ub‐path t he shortest, nei ther pro‐ v i d e t h e t o t a l p a t h l e n g t h a n d t h e b a l a n c e d e g r e e . T h e r e f o r e , w e propose t h e bal a nce rat i o (R) to measur e the workload b alance. The c a lculation formula is (longest sub ‐path length ‐ shortest sub ‐ path length)/average length of each sub ‐path × 100% , given th e t o tal path l ength of e ach stand‐ a r d d a t a s e t i n t h e c a s e o f d i f f e r e n t t r a v e l e r n u m b e r s a n d t h e calculation r esults o f balance de‐ gree are sho wn in Ta ble 5. I t c a n b e s e e n f r o m T a b l e s 3 t o 5 t h a t i n c o m p a r i s o n w i t h o t h e r a l g o r i t h m s , t h e V N S ‐ G A c a n e n s u r e t h a t t h e v a l u e o f t h e t o t a l p a t h i s w i t h i n a n a c c e p t a b l e r a n g e w h i l e e n s u r i n g t h a t t h e v a l ‐ ue o f the longest sub‐path i s sh ort e st. Thus, it c an b e s een a s V NS‐GA can fi nd s olutions w ith b e t‐ ter balance and shorte r total path length. Fig. 6 is a det a iled path di agra m of the optim al solu t ion o f t h e VNS‐GA algori t hm aft er ru nning o n e a c h d a t a s e t . I t c a n b e i n t u i t i v e l y s e e n f r o m th e fi gure t ha t t h e n u m b e r o f n o d e s c o n t a i n e d i n each p ath is r oughly e quivalent, a nd t he travel distance r equir ed f or e ach t r avelin g sale sman i s roughly the same. I n t h e c a s e o f d i f f e r e n t n u m b e r s o f t r a v e l e r s , t h e V N S ‐ G A h a s a c h i e v e d g o o d r e s u l t s w h e n e x ‐ perimenting on several standard d at a sets. It indic ates th a t th e VNS can b e tt er solve the situation w h e r e G A i s e a s y t o f a l l i n t o l o c a l o p t i m a l i t y . V N S h a s a d e e p e r search o f the solution s pace, and the insertio n operator d esigned in th e n eighbourh ood a ction c a n o pti m ize the worklo ad b alanc e problem so t hat th e s o lut i on can be d evelop ed i n the direction o f e q u i l i b r i u m a n d t h e t o t a l p a t h is shorter.  Table 5 Balanc e of VNS‐GA on standard data set Data n m Total path Longest path R 51 51 3 495 165 0 % 51 5 552 121 29.04 % 51 10 1100 112 4.5 % 100 100 3 25524 8613 2.60 % 100 5 31222 6445 9.40 % 100 10 49606 5764 22.48 % 100 20 101235 5395 26.90 % 150 150 3 31960 10878 4.76 % 150 5 35754 7711 17.98 % 150 10 55505 5937 11.03 % 150 20 102650 5754 14.50 % Fig. 6 Path diagram of VNS‐GA on standard data set Improved Genetic Algorithm (VNS-GA) using polar coordinate classification for workload balanced multiple Traveling … 5. Conclusion In current times, the research on TSP has been relatively mature, but the mTSP of workload bal- ance involves more constraints and the problem is more complicated. Thus, there are relatively few studies. This paper proposes an improved genetic algorithm (VNS-GA) to solve mTSP with workload balance. Aiming at the goal of workload balance in the mTSP, firstly, the polar coordi- nate classification algorithm is designed to reduce the intersection between paths and quickly obtain better initial solutions. Then, a distance comparison insertion operator is designed, which specifically takes out the node in the longest path and inserts it into the shortest path to achieve the goal of workload balance faster and to improve the efficiency of the algorithm. To avoid the genetic algorithm falling into the local optimal solution prematurely, the variable neighbour- hood descent process is introduced to generate offspring, and different neighbourhood actions are used to search alternately the neighbourhood solution space adequately. Finally, the algo- rithm is tested on the standard data set of TSPLIB. The experimental results showed that the improved genetic algorithm (VNS-GA) is very competitive. VNS-GA is superior to other improved genetic algorithms on small and medium-sized example sets especially when the number of travelers is small. At the same time, there are still a lot of crossed paths in the detailed path graph, which indi- cates that the algorithm still has room for improvement. In the future, other neighbourhood ac- tions with better performance can be introduced to increase the space for the search of the algo- rithm and to improve the search efficiency. Secondly, how to find a better solution when the number of travelers is large remains to be studied further. In addition, the performance of the algorithm on very large data sets still needs further study. Acknowledgement We thank the Editor and two anonymous referees for their many helpful comments on an earlier version of our paper. This work was supported in part by the National Natural Science Foundation of China under grant numbers 72171016; and Beijing Social Science Foundation under grant numbers 20JCC005; and the Beijing Logistics Informatics Research Base. References [1] Bektas, T. (2006). The multiple traveling salesman problem: An overview of formulations and solution proce- dures, Omega, Vol. 34, No. 3, 209-219, doi: 10.1016/j.omega.2004.10.004. [2] Bostanci, B., Karaağaç, A. (2019). Investigating the shortest survey route in a GNSS traverse network, Tehnički Vjesnik – Technical Gazette, Vol. 26, No. 2, 355-362, doi: 10.17559/TV-20170924174221. [3] Iqbal Ali, A., Kennington, J.L. (1986). The asymmetric M-traveling salesmen problem: A duality based branch- and-bound algorithm, Discrete Applied Mathematics, Vol. 13, No. 2-3, 259-276, doi: 10.1016/0166-218X(86) 90087-9. [4] Gavish, B., Srikanth, K. (1986). An optimal solution method for large-scale multiple traveling salesman problems, Operations Research, Vol. 34, No. 5, 698-717, doi: 10.1287/opre.34.5.698. [5] Yu, Q.S., Lin, D.M., Wang, D. (2012). An overview of multiple traveling salesman problem, Value Engineering, Vol. 31, No. 2, 166-168, doi: 10.14018/j.cnki.cn13-1085/n.2012.02.143. [6] Carter, A.E., Ragsdale, C.T. (2006). A new approach to solving the multiple traveling salesperson problem using genetic algorithms, European Journal of Operational Research, Vol. 175, No. 1, 246-257, doi: 10.1016/j.ejor.2005. 04.027. [7] Singh, A., Baghel, A.S. (2009). A new grouping genetic algorithm approach to the multiple traveling salesperson problem, Soft Computing, Vol. 13, 95-101, doi: 10.1007/s00500-008-0312-1. [8] Zhou, W., Li, Y. (2010). An improved genetic algorithm for multiple traveling salesman problem, In: Proceedings of 2 nd International Asia Conference on Informatics in Control, Automation and Robotics (CAR 2010), Wuhan, Chi- na, 493-495, doi: 10.1109/CAR.2010.5456787. [9] Koh, S.P., bin Aris, I., Ho, C.K., Bashi, S.M. (2006). Design and performance optimization of a multi-TSP (Traveling Salesman Problem) algorithm, Artificial Intelligence and Machine Learning AIML, Vol. 6, No. 3, 29-33. [10] Hu, S.J., Lu, H.Y., Huang, Y., Xu, K.B. (2019). Improved genetic algorithm for solving multiple traveling salesman problem with balanced workload, Computer Engineering and Applications, Vol. 55, No. 17, 150-155. [11] Guo, S. (2019). Solutions space analysis of MTSP and application in VRP optimization, Beijing University of Posts and Telecommunications, Beijing, China, from https://kns.cnki.net/KCMS/detail/detail.aspx?dbname=CMFD201902&filename=1019113259.nh, accessed April 11, 2021. Advances in Production Engineering & Management 16(2) 2021 183 Wang, Lu, Shen [12] Lu, Z., Zhang, K., He, J., Niu, Y. (2016), Applying K-means clustering and genetic algorithm for Solving MTSP, In: Gong, M., Pan, L., Song, T., Zhang, G. (eds.), Bio-inspired Computing – Theories and Applications, Springer Singa- pore, 278-284, doi: 10.1007/978-981-10-3614-9_34. [13] Cheikhrouhou, O., Khoufi, I. (2021). A comprehensive survey on the multiple travelling salesman problem: Appli- cations, approaches and taxonom, Computer Science Review, Vol. 40, doi: /10.1016/j.cosrev.2021.100369. [14] Pan, J., Wang, D. (2006). An ant colony optimization algorithm for multiple travelling salesman problem, In: Proceedings of the First International Conference on Innovative Computing, Information and Control - Volume I (ICICIC'06), Beijing, China, 210-213, doi: 10.1109/icicic.2006.40. [15] Venkatesh, P., Singh, A. (2015). Two metaheuristic approaches for the multiple traveling salesperson problem, Applied Soft Computing, Vol. 26, 74-89, doi: 10.1016/j.asoc.2014.09.029. [16] Ryan, J.L., Bailey, T.G., Moore, J.T., Carlton, W.B. (1998), Reactive tabu search in unmanned aerial reconnaissance simulations, In: Proceedings of the 30 th 1998 Winter Simulation Conference. Proceedings (Cat. No.98CH36274), Washington, USA, Vol. 1, 873-879, doi: 10.1109/wsc.1998.745084. [17] Song, C.H., Lee, K., Lee, W.D. (2003). Extended simulated annealing for augmented TSP and multi-salesmen TSP, In: Proceedings of the International Joint Conference on Neural Networks 2003, Oregon, USA, Vol. 3, 2340-2343, doi: 10.1109/IJCNN.2003.1223777. [18] Hu, Y., Yao, Y., Lee, W.S. (2020). A reinforcement learning approach for optimizing multiple traveling salesman problems over graphs, Knowledge-Based Systems, Vol. 204, Article No. 106244, doi: 10.1016/j.knosys.2020. 106244. [19] Bonz, J. (2021). Application of a multi-objective multi traveling salesperson problem with time windows, Public Transport, Vol. 13, 35-57, doi: 10.1007/s12469-020-00258-6. [20] Liu, H., Zhang, H., Xu, Y. (2021). The m-Steiner traveling salesman problem with online edge blockages, Journal of Combinatorial Optimization, Vol. 41, 844-860, doi: 10.1007/s10878-021-00720-6. [21] Dong, H.Y., Huang, M., Wang, X.W., Zheng, B.L. (2009). Review of variable neighborhood search algorithm, Control Engineering of China, Vol. 16, No. 2, 1-5. [22] Li, J., Guo, Y.H. (2001). Theory and method of optimal scheduling of logistics distribution vehicles, China Fortune Press, Beijing, China. [23] Min, J.N., Jin, C., Lu, L.J. (2019). Split-delivery vehicle routing problems based on a multi-restart improved sweep approach, International Journal of Simulation Modelling, Vol. 18, No. 4, 708-719, doi: 10.2507/IJSIMM18(4)CO19. [24] Min, J.N., Jin, C., Lu, L.J. (2019). Maximum-minimum distance clustering method for split-delivery vehicle-routing problem: Case studies and performance comparisons, Advances in Production Engineering & Management, Vol. 14, No. 1, 125-135, doi: 10.14743/apem2019.1.316. [25] Xiong, C., Wu, H.P., Li, B. (2010). Improved genetic algorithm for solving MTSP, In: Proceedings of the 4 th China Intelligent Computing Conference, Beijing, China, 143-149. [26] Hu, S.J. (2019). Research on multiple traveling salesman problem based on improved genetic algorithm, Master Thesis, Jiangnan University, Jiangnan, China. 184 Advances in Production Engineering & Management 16(2) 2021