% random_search.m -- Genetic algorithm, greedy algorithm, simulated annealing % and random search algorithms demonstration for solving TSP problem % N: number of cities. Cities are labeled from 1 to N % taking advantage of rotation invariance, we start every tour at city #1 % we did not take advantage of symmetricity yet % implemented algorithms: % 1. exhaustive search (only if N<= 8) % 2. genetic algorithm (using tsp_ga.m) % 3. greedy search (Lin and Kernighan) % 4. random search % 5. simulated annealing % call rshift.m, permugen.m % % (C) 2001-2010 by Yu Hen Hu % created: 12/1/2001, last modified: Nov. 29, 2010 clear all, close all, max_ncity=8; % max. # of cities that exhaustive search will be performed. N=input(['# of cities (default = ' int2str(max_ncity) ') = ']); if isempty(N), N=max_ncity; end if N<=max_ncity, disp(['Exhaustive search will be performed when # of cities <= ' ... int2str(max_ncity) '.']); else disp('Exhaustive search will NOT be performed due to too many cities!'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (1) generate a TSA problem for N cities pos=rand(N,2); % cities located within unit square d_cities=dist(pos,pos'); % matlab dist function first matrix W is row vector % second P matrix counts column vector dist(W,P) figure(1) subplot(231),plot(pos(:,1),pos(:,2),'o') axis([0 1 0 1]), axis square, hold on for i=1:N text(pos(i,1)+.03,pos(i,2),int2str(i)); end hold off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % exhaustive search disp('********** EXHAUSTIVE SEARCH METHOD'); % costmin=5; % initialize with an upper bound of cost if N <=max_ncity, tmp=permugen([2:N]); alltour=[ones(size(tmp,1),1) tmp ]; % each row is a tour for i=1:factorial(N-1), cost(i)=sum(diag(d_cities(alltour(i,:)',rshift(alltour(i,:))'))); end [costmin,idx]=min(cost); tourmin=alltour(idx,:); disp(['cost function evaluation: ' int2str(factorial(N-1)) ' times!']) disp(['minmum trip length = ' num2str(costmin)]) disp('optimum tour = ') disp(num2str(tourmin)) subplot(232),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on trip=pos(tourmin',:); trip=[trip; trip(1,:)]; plot(trip(:,1),trip(:,2),'-') title(['Exhaustive search, cost = ' num2str(costmin)]) hold off end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Genetic search % a public domain matlab mfile tsa_ga.m is used. % tsa_ga.m: author by Joseph Kirk, posted on 16 Jan 2007 (Updated 02 Jun 2009) % http://www.mathworks.com/matlabcentral/fileexchange/13680-traveling-sales % man-problem-genetic-algorithm % Input: % XY (pos) is an Nx2 matrix of city locations, where N is the number % of cities % DMAT (d_cities) is an NxN matrix of point to point distances/costs % POP_SIZE (ngpool) is the size of the population (should be divisible by 4) % NUM_ITER (ngen) is the number of desired iterations for the algorithm to run % SHOW_PROG (default=0) shows the GA progress if true % SHOW_RES (default=0) shows the GA results if true % disp('********** GENETIC ALGORITHM (tsa_ga.m)'); ngen=input('# of generations to evolve (default = 10000) = '); if isempty(ngen), ngen=1e4; end ngpool=input('# of chromosoms in the gene pool (default = 60) = '); % size of gene pool if isempty(ngpool), ngpool=60; end % gpool=zeros(ngpool,N); % gene pool % for i=1:ngpool, % intialize gene pool % gpool(i,:)=[1 randomize([2:N]')']; % end show_prog = 0; show_res = 0; [tourmin,costmin] = tsp_ga(pos,d_cities,ngpool,ngen,show_prog,show_res); % costmin=N, tourmin=zeros(1,N); cost=zeros(1,ngpool); % for i=1:ngen, % % step 1. evaluate fitness of current chromosoms % % for TSP problem, it is the trip length. The small the better % for i=1:ngpool, % cost(i)=sum(diag(d_cities(gpool(i,:)',rshift(gpool(i,:))'))); % end % % record current best solution % [costmin,idx]=min(cost); % tourmin=gpool(idx,:); % % % step 2. cross breeding and mutation % % since each chromosom is an integer vector, cross-breeding can not be % % done by simply combining binary vectors. % % Instead, we let the off-spring to keep the common genes of parents, and % % randomly shuffle genes when parents disagree % % for simplicity, we only cross-breed the two best solutions % % [csort,ridx]=sort(cost); % sort from small to big. % father=gpool(ridx(1),:); mother=gpool(ridx(2),:); % parents % % find index of genes that are the same in father and mother % sameidx=[father==mother]; % include index #1, a mask vector % diffidx=find(sameidx==0); % idx of different elements % if length(diffidx) <=4, % father and mother too close! % % do mutation % boy=[1 randomize([2:N]')']; % girl=[1 randomize([2:N]')']; % else % boy=father.*sameidx; boy(diffidx)=randomize(father(diffidx)')'; % girl=mother.*sameidx; girl(diffidx)=randomize(mother(diffidx)')'; % end % gpool=[gpool(ridx(1:ngpool-2),:);boy;girl]; % replace the worst two % end disp(['cost function evaluation: ' int2str(ngen*2+ngpool) ' times!']) disp(['minmum trip length = ' num2str(costmin)]) disp('optimum tour = ') disp(num2str(tourmin)) subplot(233),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on trip=pos(tourmin',:); trip=[trip; trip(1,:)]; plot(trip(:,1),trip(:,2),'-') title(['Genetic Search, cost = ' num2str(costmin)]) hold off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lin and Kernighan hill climbing result % Use local interchange to search for steepest descending local minimum solution % disp('********** GREEDY SEARCH METHOD'); nn=input('local neighborhood size (default = 2) = '); if isempty(nn), nn=2; end tour0=[1 randomize([2:N]')']; % random initial tour locmin=N; converged = 0; iter=0; while converged == 0, iter=iter+1; % generate neighborhood trips tmp=nbgen(tour0(2:N),nn); nbtour=[tour0; [ones(size(tmp,1),1) tmp]]; % all the neighboring tours and initial tour ntour=size(nbtour,1); cost=zeros(1,ntour); for i=1:ntour, cost(i)=sum(diag(d_cities(nbtour(i,:)',rshift(nbtour(i,:))'))); end [locost,idx]=min(cost); % 1 step steepest descent search result tour=nbtour(idx,:); if locost < locmin, % if still improvement, continue to next iteration locmin=locost; loctour=tour; tour0=tour; elseif locost == locmin, % a local minimum is reached converged=1; costmin=locost; tourmin=tour; end end ncomplexity=iter*prod([N-1:-1:N-nn])/factorial(nn); disp(['cost function evaluation: ' ... int2str(ncomplexity) ' times!']) disp(['minmum trip length = ' num2str(costmin)]) disp('optimum tour = ') disp(num2str(tourmin)) subplot(234),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on trip=pos(tourmin',:); trip=[trip; trip(1,:)]; plot(trip(:,1),trip(:,2),'-') title(['Greedy search, cost = ' num2str(costmin) ', nb size = ' int2str(nn)]) hold off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % randomized search disp('********** RANDOM SEARCH METHOD'); ntrial = input('max # of search steps (default=10000) = '); if isempty(ntrial), ntrial=1e4; end tour0=[1 randomize([2:N]')']; % initial tour costmin=N; for nt=1:ntrial, cost(nt)=sum(diag(d_cities(tour0',rshift(tour0)'))); if cost(nt) costmin, % if no reduction % toss a dice dice=rand; % range between 0 and 1 if dice <= T(nt), % accept higher cost tour0=tour; end end % convergence test if nt == ntrial | ccount== ccountmax, converged=1; end end disp(['cost function evaluation: ' int2str(nt) ' times!']) disp(['minmum trip length = ' num2str(costmin)]) disp('optimum tour = ') disp(num2str(tourmin)) subplot(236),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on trip=pos(tourmin',:); trip=[trip; trip(1,:)]; plot(trip(:,1),trip(:,2),'-') title(['SA search, cost = ' num2str(costmin)]) hold off