Genetic Algorithm for VRP

来源：互联网发布：js 定义json字符串编辑：程序博客网时间：2024/05/18 02:28

Part 1: Representation and explanation of the given problem

This is a classical optimization problem in graphic thoery. Variable descriptions and mathmatical model are given below:

Variables:

K: the number of vehicles;

Q: capacity of each vehichle (all vehicles have the same capacity);

k: index of each vehicle;

S_k: length travelled by the kth vehicle;

N_k: the number of customers assigned to the kth vehicle;

M: the number of all customers;

x_mk: x coordination of the mth customer corresponding to the kth vehicle;

y_mk: y coordination of the mth customer corresponding to the kth vehicle;

D_mk: demand amount of the mth customer corresponding to the kth vehicle;

T_mk: arriving time of the mth customer corresponding to the kth vehicle;

T¹_mk: ready time of the mth customer corresponding to the kth vehicle;

T²_mk: due time of the mth customer corresponding to the kth vehicle;

T_k: the time when the kth vehicle returns to the original starting location;

T²₀: the total consuming time limit;

t: duration of sevice time ( same in all locations );

S: total distance of all vehicles;

n: the orignal index of the given location.

d_k(i , j): distance between the ith and jth location for the kth vehicle;

Analysis:

The goal is to fulfill demands of all customers at the lowest sum of travelling distances by all vehicles according to the problem requirement. For operational sake, all the customers have to be grouped into K sections. The kth section is supposed to contain N_k specific customers which are assigned to the kth vehicle.

Constraints can be listed in two aspects:

1. number and volumes of the vehicles;

2. time constraints including: ready time and due date of each customer and the total amount of time used by each vehicle. The latter one is mentioned in files as due date of the 0^thlocation.

Mathematical model:

Objective:

in which, we give the general mathematical notation of the optimization objective and overall restraints. Some detailed relations among these varialbes are stated as following:

where d(m,m+1) originally represents the distance between the mth and (m+1)th customer, yet the task description thinks that this is eaqual to the time consumed in the jounry between the two locations. As a result, we directly adopt d(m,m+1) as the time. t is the service time in the mth location for kth vehicle. The total time during the jouney of the kth vehicle T_k can be seen as the sum of 3 parts: 1 elasped time when the kth vehicle arrives the last customer N_k; 2 sevice duration at the N_kth location; 3 time used in the way from the N_kth location to the starting point.

Part 2: Genetic Algorithm for VRP

The mathematical model only lists the general relations among all the variables and give a overall description of the problem. However, several critical issues call for solutions:

1. How to group these customers?

2. For a certain group, how to choose the exact route?

In fact, the two issues are mutually reliant. If the former one can be solved, then VRP problem can be transformed into the TSP (Travelling Salesman Problem). Distances in certain routes will definitely influence the way to devide customers.

Chromosome generating:

The coding of VRP problem is quite different from norm problems. All the original chromosomes must obey all the constraints above, as a result, we are confromted with some chromosome structure problem.

1. Time restraint of all customers. We build a directional graph spanned by all the locations. We know that the p^th customer has a ready time T¹_p and a due time T²_p, there is a route from customer p to customer q only when

which describs that if there is a vehicle moves from customer p, with time consumed on the way to q^th customer and service time at location p taken into consideration, the time when the vehicle arrives at the q^th location must lie between the ready time and due time of the the q^th customer.

By this, we successfully transfer a complete graph into a directional graph whose degree is substantially decreased, and many routes need not to be considered.

2. Time constraint of a route. As is stated above, the time of traversing all the locations of a route and returning to the starting point must be less than the due time of the starting point. By further examination of some given examples and routine thinking we know that if there is no violation of the time restraint of the last customer of a route, the vehicle can return to the starting point on time.

3. Capacity limitation. In detailed operation, we will randomly choose a customer as the first customer of a certain route, and then there only some adjacent points next to it, from which we choose one as the second customer of the route, and so on. At each location, we will compute the cumulative demand of the customers on present route, a cease of the route when the cumulative demand exceeds the vehicle capacity. Then choose another customer as the first point of a new route and excute the same procedure. It must be addressed that each customer can only be contained in one route.

4. Constaint of the number of vehicles. We can excute the progress stated in Step 3, if the chromosome uses more vehicles than allowed, regenerate a new one until it satisfy all vehicle number constraint.

Estimation function:

Suppose there are M customers and K vehicles, there should be M+K+1 elements in the serial above. Our objective is to minimize the total distance, which in the specific serial can be expressed as the sum of all distances between two consecutive elements. I would like to choose such a calculation as fitness value.

where i is the index in the serial and d(i, i+1) represents the distance between the i^thand (i+1)^th location.

Crossover and variation:

Crossover here are also barried by the structual problem (all the constraints), we can not merely choose a crossover point and exchange haves of two chromosomes as normal. As the problem is based on the directional graph, and then we need only think over constraints of capacity limitation and vehicle number.

We randomly choose two chromosomes from parental generation. We build a new directional graph among customers only(no consideration of the starting point), all the routes between points are selected from the two chromosomes. In fact, quite a lot are repeated in two chromosomes, each vertex has no more than two inputs and two outputs. Maybe some points are isolated (in two chromosomes the point is a route). Then we can simplify the directional graph. We firstly choose a shorter input of two-inputs point and desert the other one. Then choose a shorter output of two-outputs and let the other abandoned. With the redued graph, we furtherly generate the vehicle routes as above, in which capacity are taken account.

In each ieration, we randomly generate several chromosomes as variation.

Part3: Calculation result:

We choose the data listed in file “C1_2_1.txt” as an example.

Each generation 10 chromosomes are selected and 3 iterations are taken. We only give the exact chromosome with the best fitness.

Initial generation:

Index

Fitness

(min)

4292.9

4601.7

5508.7

4791.8

4391.6

7438.8

7361.2

6563.2

5972.3

6705.4

Best fitness: 4292.9

Best chromosome:

Index

Route

0 20 41 12 129 11 6 122 90 67 17 39 0

0 21 23 182 75 163 194 145 195 52 92 0

0 24 61 139 0

0 30 120 19 192 196 97 14 89 105 15 59 0

0 32 171 65 86 115 94 51 110 162 0

0 45 155 0

0 144 119 166 35 126 71 9 1 99 53 0

0 100 64 179 109 149 56 0

0 57 5 10 193 46 128 106 167 34 95 158 0

0 114 159 38 150 22 151 16 140 187 142 111 63 0

0 73 116 147 160 47 91 70 0

0 78 175 13 43 2 0

0 190 82 0

0 85 80 31 25 172 77 0

0 93 55 18 54 185 132 7 181 188 108 0

0 101 0

0 133 48 26 152 40 153 169 96 130 28 74 107 0

0 118 83 143 176 36 33 121 165 49 0

0 148 103 197 124 141 69 200 174 136 189 0

0 161 104 0

0 164 66 0

0 170 134 50 156 112 168 79 29 87 42 123 198 0

0 177 3 88 8 186 127 98 157 0

0 178 27 173 154 0

0 135 58 184 199 37 81 138 137 183 0

0 180 84 191 125 4 72 0

0 44 102 146 0

0 60 0

0 62 0

0 68 76 0

0 113 0

0 117 0

0 131 0

Second generation:

Index

Fitness

(min)

3559.5

3410.8

3712.5

4246.4

3852.1

6075.3

6492.6

7536.1

6920.8

5671.4

Best fitness: 3410.8

Best chromosome:

Index

Route

0 20 41 85 80 31 25 172 77 110 162 0

0 21 23 182 75 163 194 145 195 52 92 0

0 30 120 19 192 196 97 14 96 130 28 74 149 0

0 32 171 65 86 115 94 51 174 136 189 0

0 45 27 173 154 24 61 100 64 179 109 0

0 161 104 18 54 185 132 7 181 117 49 0

0 101 144 119 166 35 126 71 9 1 99 53 0

0 50 156 112 168 79 29 87 42 123 198 56 0

0 57 118 83 143 176 36 33 121 165 188 108 0

0 133 48 26 152 40 153 169 89 105 15 59 0

0 60 82 180 84 191 125 4 72 0

0 62 131 44 102 146 68 76 0

0 38 150 22 151 16 140 187 142 111 63 0

0 93 55 135 0

0 113 155 78 175 13 43 2 90 67 17 39 107 0

0 116 12 129 11 6 122 139 0

0 148 103 197 124 141 69 200 0

0 190 5 10 193 46 128 106 167 34 95 158 0

0 164 66 147 160 47 91 70 0

0 170 134 0

0 177 3 88 8 186 127 98 157 0

0 178 0

0 58 184 199 37 81 138 137 183 0

0 73 0

0 114 0

0 159 0

Third generation:

Index

Fitness

(min)

2877.3

3602.6

3166.5

4982.7

4721.6

6282.7

7084.9

6716

6402.1

6907.6

Best fitness: 2877.3

Best chromosome:

Index

Route

0 20 41 85 80 31 25 172 77 110 162 0

0 21 23 182 75 163 194 145 195 52 92 0

0 30 120 19 192 196 97 14 89 105 15 59 0

0 32 171 65 86 115 94 51 174 136 189 0

0 45 27 173 154 24 61 100 64 179 109 149 0

0 161 104 18 54 185 132 7 181 117 49 0

0 101 144 119 166 35 126 71 9 1 99 53 0

0 93 55 135 58 184 199 37 81 138 137 183 56 0

0 57 118 83 143 176 36 33 121 165 188 108 0

0 60 82 180 84 191 125 4 72 0

0 62 131 44 102 146 68 76 96 130 28 74 0

0 114 159 38 150 22 151 16 140 187 142 111 63 0

0 164 66 147 160 47 91 70 0

0 73 116 12 129 11 6 122 139 0

0 113 155 78 175 13 43 2 90 67 17 39 107 0

0 133 48 26 152 40 153 169 0

0 148 103 197 124 141 69 200 0

0 190 5 10 193 46 128 106 167 34 95 158 0

0 170 134 50 156 112 168 79 29 87 42 123 198 0

0 177 3 88 8 186 127 98 157 0

0 178 0

The example given above shows that performance improved rapidly with the iteration of genetic algorithm. By perceptual knowledge, we manage it by two ways: (1) In crossover progress, we cut off the longer edge of the two-input or two-output vertex to keep it connected to other vertexes by the remained one. We see that the total objective function if the distance of all vehicle routes, when each element of the vehicle route becomes shorter, there must be an reduction in the total distace. (2) Each vehicle route has two extented cost of distance: the first customer and the last customer of the vehicle route connected to the starting point. As the iteration going on, the number of vehicles really used decreases, so the distance reduction in the extended cost will ensure the total distance optimized.

This is a strictly constrainted problem, which makes it hard to variate. We directly randomly select some new chromosomes, in which some new genes are included. It has reached the rule of traditional rule of getting new genes. We think that it makes sense.

The following gives the matlab codes:

You only need to excute :

[resultChrom,lastMin,fitAll,ddd]=geneticVRP(iter)

in the command window, where iter denotes the iteration of genetic method.

% geneticVRP.m

function [resultChrom,lastMin,fitAll,ddd]=geneticVRP(iter)

clear

iter=4;

% Assumptions:

% (1) service time of all customers is the same

% (2) time and demand constraints can never be violated

% (3) K as the number of vehicles is only a upper limit, we need not use

% all the vehicles

% (4) no location but the start point can be revisited.

%55555555555555555555555555555555555555555555555555555555555555555555555555

% Some preparations

%55555555555555555555555555555555555555555555555555555555555555555555555555

%5 parameter description

originData11; % load data

[height,width]=size(originData);

M=height-1; % number of customers

K=50; % number of vehicles

capacity=200; % capacity of each vehicle

NC=5;

coord=originData(2:M+1,2:3); % coordinations

timLim=originData(2:M+1,5:6); % time limits

demand=originData(2:M+1,4); % demand

serTime=originData(2,7); % service time

%5 calculate distances of all points to start point

disOrg=zeros(1,M);

for i=1:M

disOrg(i)=sqrt((originData(i+1,2)-originData(1,2))^2+(originData(i+1,3)-originData(1,3))^2);

end

clear originData;

%5 build directional graph

pt2={};

ctPt=zeros(1,M);

dis2=ones(M)*Inf;

for i=1:M-1

for j=i+1:M

disBetween=sqrt((coord(i,1)-coord(j,1))^2+(coord(i,2)-coord(j,2))^2)+serTime; % distance between two consecutive locations

if(timLim(i,2)+disBetween>timLim(j,1)&&timLim(i,1)+disBetween<timLim(j,2))

ctPt(i)=ctPt(i)+1;

pt2{i}(ctPt(i))=j;

dis2(i,j)=disBetween;

end

if(timLim(j,2)+disBetween>timLim(i,1)&&timLim(j,1)+disBetween<timLim(i,2))

ctPt(j)=ctPt(j)+1;

pt2{j}(ctPt(j))=i;

dis2(j,i)=disBetween;

end

%55555555555555555555555555555555555555555555555555555555555555555555555555

% Genetic job

%55555555555555555555555555555555555555555555555555555555555555555555555555

% choose nc couples of chromosomes which fulfill all the constraints

% calculate the fitness of all the chromosomes

chrom={};

numVehi=zeros(1,NC*2);

fit=ones(1,NC*2);

for i=1:NC*2

while(1)

[chrom{i},numCus]=randomSelect(M,K,demand,disOrg,ctPt,pt2,dis2,capacity);

if(numCus==M)

break;

end

% chromosome crossover

% let the NC couples of chromosomes randomly unioned

% regenerate some chromosomes and then select the best 2*NC chromosomes as

% the parental chromosomes

% Get a discrete seria

chromChild={};

ddd={};

fitAll=zeros(iter,2*NC);

lastMin=inf;

lastIndex=0;

for j=1:iter

index=randperm(2*NC);

for i=1:NC

chromChild{i}=crossOver(chrom{index(2*i-1)},chrom{index(2*i)},dis2,M,K,capacity,demand);

end

for i=1:NC

chrom{i}=chromChild{i};

while(1)

[chrom{i+NC},numCus]=randomSelect(M,K,demand,disOrg,ctPt,pt2,dis2,capacity);

if(numCus==M)

break;

end

fit(i)=fitness(chrom{i},dis2-serTime,disOrg);

fit(i+NC)=fitness(chrom{i+NC},dis2-serTime,disOrg);

end

[minChrom,minIndex]=min(fit);

if(minChrom==lastMin)

break;

end

lastMin=minChrom;

lastIndex=minIndex;

lend=length(chrom{lastIndex});

ddd{j}=[];

for k=1:lend

ddd{j}=[ddd{j},0,chrom{lastIndex}{k}];

end

fitAll(j,:)=fit;

end

resultChrom=chrom{lastIndex};

%originalData11.m

originData=[0 70 70 0 0 1351 0

1 33 78 20 750 809 90

2 59 52 20 553 602 90

3 10 137 30 147 219 90

4 4 28 10 616 661 90

5 25 26 20 128 179 90

6 86 37 10 478 531 90

7 1 109 10 616 680 90

8 6 135 40 351 386 90

9 32 79 20 655 721 90

10 24 26 20 219 271 90

11 86 36 20 384 443 90

12 95 35 10 196 252 90

13 63 50 10 364 416 90

14 100 106 10 567 626 90

15 99 112 20 846 911 90

16 36 135 10 598 669 90

17 57 59 10 824 890 90

18 8 124 10 231 298 90

19 85 106 20 186 253 90

20 103 69 30 33 93 90

21 109 131 20 72 141 90

22 43 140 10 410 470 90

23 115 134 30 140 199 90

24 98 70 10 467 511 90

25 112 67 10 462 534 90

26 102 104 20 201 258 90

27 93 75 30 178 238 90

28 90 104 20 845 910 90

29 127 108 10 677 743 90

30 84 99 20 32 102 90

31 113 69 20 370 440 90

32 129 9 10 84 147 90

33 18 38 30 582 638 90

34 30 27 10 772 835 90

35 25 80 20 380 440 90

36 17 37 30 493 544 90

37 32 106 10 602 657 90

38 43 135 10 215 293 90

39 61 59 20 931 971 90

40 104 106 20 390 434 90

41 109 71 20 95 163 90

42 121 110 30 870 923 90

43 61 48 20 448 517 90

44 74 99 20 176 243 90

45 89 73 10 19 94 90

46 21 25 20 405 453 90

47 99 28 30 392 442 90

48 101 96 10 101 160 90

49 9 114 20 904 954 90

50 121 112 20 221 277 90

51 137 6 10 611 678 90

52 118 131 20 782 857 90

53 34 82 10 928 1000 90

54 4 125 30 319 398 90

55 26 105 40 115 179 90

56 35 123 20 1130 1198 90

57 21 48 10 53 123 90

58 21 115 10 306 376 90

59 99 108 20 952 992 90

60 1 34 20 77 132 90

61 94 70 20 556 610 90

62 74 93 20 23 66 90

63 32 128 20 1079 1123 90

64 94 73 10 741 791 90

65 135 3 20 250 300 90

66 97 28 20 110 170 90

67 59 56 20 745 780 90

68 80 97 10 457 518 90

69 134 31 20 508 566 90

70 103 34 30 577 636 90

71 30 79 20 561 632 90

72 4 30 40 693 768 90

73 87 39 20 35 92 90

74 90 101 10 933 1008 90

75 119 136 30 334 384 90

76 78 92 30 555 610 90

77 110 66 10 650 710 90

78 62 58 20 173 229 90

79 127 112 20 595 636 90

80 113 71 10 286 341 90

81 34 104 10 683 760 90

82 2 28 10 148 199 90

83 17 43 10 205 279 90

84 1 24 10 333 385 90

85 112 71 20 193 252 90

86 135 0 10 331 405 90

87 125 108 20 778 825 90

88 8 138 30 234 316 90

89 103 111 20 668 719 90

90 57 53 10 641 697 90

91 98 31 20 480 541 90

92 116 132 10 870 953 90

93 28 99 10 51 93 90

94 137 3 20 522 581 90

95 29 28 10 865 925 90

96 94 105 10 656 728 90

97 92 107 30 465 532 90

98 10 134 10 629 677 90

99 34 79 10 841 902 90

100 94 72 20 647 703 90

101 35 76 10 35 91 90

102 76 99 20 273 330 90

103 127 28 30 128 195 90

104 8 119 20 142 197 90

105 102 110 10 754 814 90

106 24 24 10 583 647 90

107 64 61 20 1011 1077 90

108 22 38 20 950 1003 90

109 88 72 30 925 979 90

110 107 63 30 748 800 90

111 33 131 30 979 1037 90

112 123 112 20 405 459 90

113 67 55 10 15 62 90

114 39 128 10 65 133 90

115 136 0 30 426 491 90

116 93 36 20 98 166 90

117 7 115 10 801 872 90

118 21 44 20 115 180 90

119 28 74 20 188 260 90

120 86 102 10 93 159 90

121 19 37 20 665 738 90

122 84 36 10 568 626 90

123 118 108 10 972 1008 90

124 132 24 10 311 384 90

125 6 25 10 516 574 90

126 28 78 10 479 529 90

127 2 132 10 526 584 90

128 22 22 10 494 550 90

129 88 35 20 285 358 90

130 94 104 10 754 813 90

131 74 96 30 80 153 90

132 0 117 30 515 585 90

133 99 96 10 38 93 90

134 119 114 20 128 185 90

135 28 107 20 206 275 90

136 133 12 10 797 874 90

137 35 102 10 877 938 90

138 36 105 20 776 852 90

139 84 39 20 659 720 90

140 34 138 20 704 749 90

141 136 29 10 406 482 90

142 33 133 20 890 941 90

143 16 39 20 314 358 90

144 32 77 20 100 157 90

145 122 131 20 608 661 90

146 77 100 20 356 429 90

147 96 26 10 202 262 90

148 127 29 20 70 141 90

149 91 95 20 1036 1097 90

150 42 136 10 310 381 90

151 36 132 10 509 571 90

152 103 106 30 296 347 90

153 107 105 10 472 539 90

154 96 74 10 351 437 90

155 63 57 10 79 140 90

156 122 113 20 307 374 90

157 14 131 10 712 784 90

158 30 31 20 960 1016 90

159 40 132 20 137 182 90

160 99 26 40 287 363 90

161 10 118 10 76 144 90

162 102 60 30 840 900 90

163 120 136 10 417 482 90

164 97 29 30 49 119 90

165 19 39 10 764 823 90

166 27 77 10 284 349 90

167 30 24 20 673 748 90

168 126 112 10 507 543 90

169 105 110 20 567 634 90

170 118 114 10 65 126 90

171 135 6 30 153 210 90

172 111 66 20 561 617 90

173 95 77 10 278 325 90

174 139 10 30 723 756 90

175 61 54 10 262 329 90

176 16 38 20 387 466 90

177 16 138 10 86 142 90

178 91 78 10 88 141 90

179 91 71 20 829 890 90

180 0 27 30 241 292 90

181 6 113 20 720 768 90

182 119 139 20 234 298 90

183 35 99 20 976 1025 90

184 31 112 20 401 481 90

185 5 119 20 425 484 90

186 3 135 20 433 491 90

187 30 136 10 792 850 90

188 20 39 10 855 913 90

189 131 11 10 895 960 90

190 28 33 20 55 128 90

191 4 23 10 414 491 90

192 87 108 20 275 350 90

193 21 26 20 311 365 90

194 121 134 20 517 567 90

195 119 130 20 692 765 90

196 91 108 10 384 430 90

197 130 27 30 229 280 90

198 101 107 10 1038 1090 90

199 34 108 10 506 566 90

200 131 31 20 604 657 90];

% radomSelect.m

function [serial_non,ctAll]=randomSelect(M,K,demand,disOrg,ctPt,pt2,dis2,capacity)

% randomSelect is to choose an original route

% denotions of all inputs are consistent with ones in m-file of graphBuild

%55555555555555555555555555555555555555555555555555555555555555555555555555

%5 randomly generate routes

%55555555555555555555555555555555555555555555555555555555555555555555555555

flagEmp=ones(1,M);

% Get a discrete serial

dd=randperm(M);

flagRand=1;

isSuc=1;

serial_non={};

capVe=zeros(1,K);

ctAll=0; % put down customers arranged

for i=1:K

pt=dd(flagRand); % Select the first customer of a route

flagRand=flagRand+1;

while(~flagEmp(pt)&flagRand<=M)

pt=dd(flagRand);

flagRand=flagRand+1;

end

if(flagRand>M) % all customers has been arranged while some vehicles are empty

break;%%'***************************'

end

timeCur=disOrg(pt); % initialization of each vehicle

demandCur=demand(pt);

count=1; % count the vehicles used

serial_non{i}(count)=pt;

capVe(i)=demand(pt);

flagEmp(pt)=0;

ctAll=ctAll+1;

%------- ------- --------- -------- -------- -------- ------- --------

judgeFlag=1;

while(judgeFlag)

while(ctPt(pt)>0) % find an unused location "next"

%index=ceil(rand(1)*ctPt(pt));

[dsd,index]=min(dis2(pt,pt2{pt}));

if(~flagEmp(pt2{pt}(index)))

pt2{pt}(index)=pt2{pt}(ctPt(pt));

ctPt(pt)=ctPt(pt)-1;

continue;

else

break;

end

isNext=ctPt(pt);

if(isNext) % verify the current point "next"

next=pt2{pt}(index);

timeCur=timeCur+dis2(pt,next);

demandCur=demandCur+demand(next);

inFlag=demandCur<=capacity;

if(inFlag)

count=count+1;

capVe(i)=demandCur;

serial_non{i}(count)=next;

flagEmp(next)=0;

ctAll=ctAll+1;

pt=next;

else

break;

end

else

break;

end

% crossOver.m

function serial_non=crossOver(parent1,parent2,dis2,M,K,capacity,demand)

% record all the routes used in the two chromosomes

inIndex=zeros(2,M);

outIndex=zeros(2,M);

len=length(parent1);

for i=1:len

subLen=length(parent1{i})-1;

for j=1:subLen

inIndex(1,parent1{i}(j+1))=parent1{i}(j);

outIndex(1,parent1{i}(j))=parent1{i}(j+1);

end

len=length(parent2);

for i=1:len

subLen=length(parent2{i})-1;

for j=1:subLen

inIndex(2,parent2{i}(j+1))=parent2{i}(j);

outIndex(2,parent2{i}(j))=parent2{i}(j+1);

end

inSum=sum(inIndex>0);

outSum=sum(outIndex>0);

for i=1:M

if(outSum(i)==2) % location is used for pt with no ending connectin to st

if(dis2(i,outIndex(1,i))>dis2(i,outIndex(2,i)))

outIndex(1,i)=0;

else

outIndex(2,i)=0;

end

if(inSum(i)==2) % location is used for pt with no start connectin to st

if(dis2(inIndex(1,i),i)>dis2(inIndex(2,i),i))

inIndex(1,i)=0;

else

inIndex(2,i)=0;

end

% firstly find sth from one-side point, and find routes with such pt as one

% end

flagUsed=zeros(1,M);

ctAll=0;

ctVe=0;

outFlag=sum(outIndex>0);

inFlag=sum(inIndex>0);

outIndex=sum(outIndex);

inIndex=sum(inIndex);

serial_non={};

for i=1:M

if(inFlag(i)==0&~flagUsed(i)&outFlag(i))

ctVe=ctVe+1;

count=1;

serial_non{ctVe}(count)=i;

capVe(ctVe)=demand(i);

pt=i;

flagUsed(i)=1;

ctAll=ctAll+1;

while(outFlag(pt)&~flagUsed(outIndex(pt)))

cap=capVe(ctVe)+demand(outIndex(pt));

if(cap>capacity)

break;

end

capVe(ctVe)=cap;

ctAll=ctAll+1;

count=count+1;

serial_non{ctVe}(count)=outIndex(pt);

pt=outIndex(pt);

flagUsed(pt)=1;

end

elseif(outFlag(i)==0&~flagUsed(i)&inFlag(i))

ctVe=ctVe+1;

count=1;

serial_non{ctVe}(count)=i;

capVe(ctVe)=demand(i);

pt=i;

flagUsed(i)=1;

ctAll=ctAll+1;

while(inFlag(pt)&~flagUsed(inIndex(pt)))

cap=capVe(ctVe)+demand(inIndex(pt));

if(cap>capacity)

break;

end

capVe(ctVe)=cap;

ctAll=ctAll+1;

count=count+1;

serial_non{ctVe}(count)=inIndex(pt);

pt=inIndex(pt);

flagUsed(pt)=1;

end

trans=serial_non{ctVe};

for i=1:count

serial_non{ctVe}(i)=trans(count+1-i);

end

% secondly try to trsemble all the points as long as possible

% we can cut it to pieces if it disobeys the capacity constraint

for i=1:M

if(flagUsed(i))

continue;

end

pre=[];

after=i;

ptPre=i;

ptAfter=i;

while(inFlag(ptPre)&~flagUsed(inIndex(ptPre)))

pre=[pre,inIndex(ptPre)];

ptPre=inIndex(ptPre);

end

while(outFlag(ptPre)&~flagUsed(outIndex(ptAfter)))

after=[after,outIndex(ptAfter)];

ptAfter=outIndex(ptAfter);

end

len1=length(pre);

len2=length(after);

ctVe=ctVe+1;

count=1;

capVe(ctVe)=0;

for j=len1:-1:1

cap=capVe(ctVe)+demand(pre(j));

if(cap>capacity)

break;

end

serial_non{ctVe}(count)=pre(j);

capVe(ctVe)=cap;

ctAll=ctAll+1;

count=count+1;

flagUsed(pre(j))=1;

end

for j=1:len2

cap=capVe(ctVe)+demand(after(j));

if(cap>capacity)

break;

end

serial_non{ctVe}(count)=after(j);

capVe(ctVe)=cap;

ctAll=ctAll+1;

count=count+1;

flagUsed(after(j))=1;

end

% fitness.m

function fit=fitness(serial_non,distance,disOrg)

% all the constraints including time and demand have been examned

% As a result, here we only figure the total distances traveled by all

% vehicles, which is the value of our objective function

len=length(serial_non);

fit=0;

for i=1:len

fit=fit+disOrg(serial_non{i}(1))+disOrg(serial_non{i}(end));

subLen=length(serial_non{i})-1;

for j=1:subLen

fit=fit+distance(serial_non{i}(j),serial_non{i}(j+1));

end