C ********************************************************* C * This program solves TSPs using simulated annealing. * C * The annealing schedule is based on Connolly's GPSIMAN.* C * Very similar to QAP code. * C * This uses the inverse operator instead of swap * C ********************************************************* PROGRAM sa INCLUDE 'SAPARAM.INC' INTEGER terminate_on_value,terminate_on_trial PARAMETER (terminate_on_value=1,terminate_on_trial=2) REAL initial_temperature,final_temperature,beta,temp1,temp2 &,min_temperature,delmin,delmax INTEGER old_fitness,qap_order,i,fitness,seed,max_fails, &next,location1,location2,trials,num_trials,fitness1,inverse, &terminate_value,best,num_passes,num_steps,kount,num_fails, &trial_best,terminate_choice,partial_calculate_fitness,f4, &solution(max_cols),old_solution(max_cols),trial,temp,j, &best_solution(max_cols),best_trial_solution(max_cols), &num_pairs,move_list(max_pairs,2),num_cities,city_matrix( &max_cities,max_cities),c2,apos,aneg,other,orig_fitness REAL temperature,random,probability_of_acceptance,temp_best DOUBLE PRECISION RN55 LOGICAL terminate/.FALSE./ CHARACTER *30 file_name c2=0 apos=0 aneg=0 other=0 CALL get_user_parameters(file_name,num_passes,seed, &terminate_value,num_trials,terminate_choice) CALL read_tsp_file(num_cities,city_matrix,file_name) c DO i=1,num_cities c PRINT *,(city_matrix(i,j),j=1,num_cities) c ENDDO qap_order=num_cities CALL form_initial_soloution(solution,qap_order,seed) c DO i=1,qap_order c PRINT *,solution(i) c ENDDO c PRINT *,"Solution",(solution(i),i=1,qap_order) CALL copy_solution(solution,best_solution,qap_order) c PRINT *,"Past it" CALL calculate_fitness(solution,city_matrix,qap_order, &old_fitness) orig_fitness=old_fitness PRINT *,"OLD",old_fitness c CALL calculate_fitness(solution,city_matrix,qap_order, c &fitness) c PRINT *,"OLD",fitness c STOP c old_fitness=bignum fitness=old_fitness max_fails=2*qap_order**2 best=fitness trial_best=fitness initial_temperature=small_number final_temperature=small_number c min_temperature=initial_temperature num_pairs=1 CALL RNSD(seed) c PRINT *,"a1" C FORM PAIR LIST DO i=1,qap_order-1 DO j=i+1,qap_order move_list(num_pairs,1)=i move_list(num_pairs,2)=j num_pairs=num_pairs+1 c PRINT *,i,j,move_list(num_pairs,1),move_list(num_pairs,2), c & num_pairs ENDDO ENDDO c PRINT *,"a2" c HERE THE USER CAN SPECIFY TO STOP AFTER A CERTIAL NUMBER OF TRIALS C HAVE BEEN REACHED OR A CERTAIN COST HAS BEEN REACHED DO WHILE (.NOT.(terminate)) c PRINT *,"A1" CALL shuffle_list(move_list,num_pairs,seed) c PRINT *,"A2" c PRINT *,"$$$",terminate_on_value,terminate_value,fitness temp_best=initial_temperature num_steps=num_passes*qap_order c PRINT *,"+++",num_steps,num_passes,qap_order C CALCULTE TEMPERATURE DECREMENT FACTOR, BETA beta=((initial_temperature-final_temperature)/(num_steps* & initial_temperature*final_temperature)) temperature=initial_temperature num_fails=0 delmin=bignum delmax=0 trial_best=fitness PRINT *,"INIT TEMP=",initial_temperature,"FINAL TEMP =", & final_temperature next=0 old_fitness=orig_fitness fitness=orig_fitness C THIS LOOP DOES A SINGLE ANNEALING RUN FROM TINIT TO TFINAL DO kount=1,num_steps f4=fitness CALL copy_solution(solution,old_solution,qap_order) C GET NEXT PAIR TO SWAP next=next+1 IF (next.GE.num_pairs) next=1 location1=move_list(next,1) location2=move_list(next,2) DO WHILE ((location1.EQ.0).OR.(location2.EQ.0)) next=next+1 IF (next.GE.num_pairs) next=1 location1=move_list(next,1) location2=move_list(next,2) ENDDO i=max(location1,location2) DO inverse=min(location1,location2),max(location1,location2) solution(inverse)=old_solution(i) i=i-1 ENDDO IF (kount.EQ.1) THEN CALL calculate_fitness(solution,city_matrix, & qap_order,fitness) ELSE fitness=partial_calculate_fitness(fitness,location1, & location2,city_matrix,solution,old_solution,num_cities) ENDIF c CALL calculate_fitness(solution,city_matrix, c & qap_order,fitness) delmax=MAX(FLOAT(ABS(fitness-old_fitness)),delmax) IF (ABS(fitness-old_fitness).GT.small_number) delmin= & MIN(delmin,FLOAT(ABS(fitness-old_fitness))) c2=c2+1 IF (fitness.GT.old_fitness) THEN C WE NOW HAVE A WORSE SOLUTION probability_of_acceptance=RN55() C STANDARD BOLTZMAN FUNCTION IN DISGUISE IF (.NOT.(temperature*LOG(1/probability_of_acceptance) & .GT.ABS(fitness-old_fitness))) THEN num_fails=num_fails+1 other=other+1 IF (num_fails.GT.max_fails) THEN C THE NUMBER OF UPHILL REJECTIONS HAS REACHED A CRITICAL NUMBER C MUST BE UBALE TO CLIMB OUT OF LOCAL MINIMUM IF (beta.NE.0) THEN min_temperature=temperature temperature=temp_best beta=0 ENDIF num_fails=0 ELSE CALL copy_solution(old_solution,solution, & qap_order) c fitness=old_fitness fitness=f4 c CALL calculate_fitness(solution,city_matrix, c & qap_order,fitness) c IF (fitness.NE.f4) PRINT *,"here",f4,fitness ENDIF c PRINT *,"%%%",probability_of_acceptance, c & EXP(-1.0*(fitness-old_fitness)/temperature), c & fitness-old_fitness,"downhill, old=",old_fitness, c & "new=",fitness c CALL calculate_fitness(solution,city_matrix, c & qap_order,fitness) ELSE c PRINT *,"This is uphill, old=",old_fitness,"new=" c & ,fitness C WE ARE ACCEPTING THE WORSE SOLUTION num_fails=0 apos=apos+1 ENDIF ELSE c PRINT *,"Solution at",fitness,"with temp",temperature aneg=aneg+1 num_fails=0 ENDIF c IF (fitness.GT.old_fitness) num_fails=0 c num_fails=0 old_fitness=fitness temperature=temperature & /(1.0+beta*temperature) IF ((fitness.LT.trial_best).OR.(kount.EQ.1)) THEN C THIS IS THE BEST SOLUTION SO FAR FOR THIS TRIAL c PRINT *,"IMPROVED SOLUTION" c PRINT *,"Solution at",fitness,"with temp",temperature c DO i=1,qap_order c PRINT *,solution(i) c ENDDO c & ,fitness1,fitness1-fitness c & ,(solution(i),i=1,qap_order) trial_best=fitness IF ((terminate_choice.EQ.terminate_on_value).AND. & (trial_best.LE.terminate_value)) THEN C WE HAVE THE COST WE WANT SO EXIT OUT OF PROGRAM best=trial_best GOTO 92 ENDIF temp_best=temperature CALL copy_solution(solution,best_trial_solution, & qap_order) ELSE c IF (MOD(kount,qap_order).EQ.1) PRINT *,"PASS ",kount ENDIF c PRINT *,"Kount=",kount ENDDO IF ((trial_best.LE.best).OR.(trial.EQ.0)) THEN best=trial_best PRINT *,"Best so far is ",best,"at trial ",trial c CALL copy_solution(best_trial_solution,best_solution, c & num_verticies) c PRINT *,qap_order c DO i=1,qap_order c PRINT *,best_trial_solution(i) c ENDDO ENDIF CALL copy_solution(best_solution,solution,qap_order) c fitness=orig_fitness c old_fitness=orig_fitness C RECALCULATE TEMPERATURE FOR NEXT TRIAL IF (initial_temperature.LE.small_number) THEN temp1=delmin+(delmax-delmin)/10.0 c PRINT *,"<<<",temp1,delmin,delmax ELSE temp1=(initial_temperature+temp_best)/2.0 ENDIF c PRINT *,"[[[",initial_temperature,temp_best IF (final_temperature.LE.small_number) THEN temp2=delmin ELSE IF (beta.NE.0) THEN temp2=final_temperature/2.0 ELSE temp2=(final_temperature+min_temperature)/2.0 ENDIF ENDIF initial_temperature=temp1 final_temperature=temp2 trial=trial+1 92 CONTINUE IF (terminate_choice.EQ.terminate_on_value) THEN IF (best.LE.terminate_value) terminate=.TRUE. ELSE IF (trial.GE.num_trials) terminate=.TRUE. ENDIF ENDDO PRINT 93, file_name, best,num_trials,num_passes,seed 93 FORMAT ('Final ',a30,i7,i7,i7,i7) END C ************************************************************* C * This subprogram forms an initial solution generic to order* C * based problems. It randomly rearranges 1..N by perfroming* C * 2-opt swaps. * C ************************************************************* SUBROUTINE form_initial_soloution(solution,qap_order,seed) INCLUDE "SAPARAM.INC" INTEGER qap_order,seed,bit,solution(max_cols),temp,move1,move2 REAL random DO bit=1,qap_order solution(bit)=bit ENDDO DO bit=1,qap_order move1=INT(random(seed)*qap_order)+1 move2=INT(random(seed)*qap_order)+1 temp=solution(move1) solution(move1)=solution(move2) solution(move2)=temp ENDDO c PRINT *,(move_list(bit),bit=1,total_bits) RETURN END C ***************************************************************** C * The move list gives the annealing engine the thext par to swap* C * This subprogram shuffles the list randomly. Connolly found * C * this to give better performance. * C ***************************************************************** SUBROUTINE shuffle_list(move_list,num_pairs,seed) INCLUDE 'SAPARAM.INC' INTEGER move_list(max_pairs,2),num_pairs,swap,move1,move2,temp,i INTEGER seed DOUBLE PRECISION RN55,m1,m2 REAL random DO swap=1,num_pairs move1=0 move2=0 DO WHILE (move1.EQ.move2) c move1=INT(random(seed)*num_pairs)+1 c move2=INT(random(seed)*num_pairs)+1 m1=RN55() m2=RN55() move1=INT(m1*num_pairs)+1 move2=INT(m2*num_pairs)+1 ENDDO DO i=1,2 temp=move_list(move1,i) move_list(move1,i)=move_list(move2,i) move_list(move2,i)=temp ENDDO ENDDO RETURN END C ********************************************************** C * get tsp details * C ********************************************************** SUBROUTINE read_tsp_file(num_cities,city_matrix,file_name) INCLUDE "SAPARAM.INC" INTEGER city_matrix(max_cities,max_cities),num_cities,city, &city_row,city_column,city_num,i,j REAL coords(max_cities,2) CHARACTER *20 file_name OPEN (UNIT=3,FILE=file_name,STATUS='OLD') READ (3,*) num_cities DO city=1,num_cities READ (3,*) city_num,coords(city,1),coords(city,2) ENDDO DO city_row=1,num_cities DO city_column=1,num_cities city_matrix(city_row,city_column)=NINT(SQRT(ABS(coords & (city_row,1)-coords(city_column,1))**2+ABS(coords( & city_row,2)-coords(city_column,2))**2)+0.5) c IF ((city_row.NE.city_column).AND.(city_matrix c & (city_row,city_column).EQ.0)) city_matrix(city_row, c & city_column)=1 ENDDO ENDDO c DO i=1,num_cities c PRINT *,(city_matrix(i,j),j=1,num_cities) c ENDDO RETURN END C *********************************************************** C * A vector copier * C *********************************************************** SUBROUTINE copy_solution(solution,old_solution,qap_order) INCLUDE 'SAPARAM.INC' INTEGER solution(max_cols),old_solution(max_cols) INTEGER qap_order,col DO col=1,qap_order old_solution(col)=solution(col) ENDDO RETURN END C ********************************************************** C * Gets the user parameters for the SA search * C ********************************************************** SUBROUTINE get_user_parameters(file_name,num_passes,seed, &terminate_value,num_trials,terminate_choice) INCLUDE 'SAPARAM.INC' INTEGER seed,terminate_value,num_passes,num_trials, &terminate_choice CHARACTER *30 file_name PRINT *,"*******************************" PRINT *,"* Simulated Annealing for TSP *" PRINT *,"* by Marcus Randall *" PRINT *,"*******************************" PRINT * PRINT *,"What is the name of Problem file?" READ *, file_name PRINT *,"How many passes for each annealing trial?" READ *,num_passes PRINT *,"What is the seed you would like to use?" READ *,seed PRINT *,"Do you want to?" PRINT *," 1. terminate on value" PRINT *," 2. terminate on trial" READ *,terminate_choice IF (terminate_choice.EQ.1) THEN PRINT *,"What value that you want the run to stop at?" READ *,terminate_value ELSE PRINT *,"How many trials?" READ *,num_trials ENDIF RETURN END C ********************************************************* C * The old random number generator. Not very effective. * C * All calls to this function can be redirected through * C * RN55 (Knuth) * C ********************************************************* REAL FUNCTION random(seed) INTEGER seed,C1,C2,oldsed SAVE oldsed DATA oldsed /0/ PARAMETER (C1=19423,C2=811) IF (oldsed.EQ.0) oldsed=seed oldsed=MOD(C1*oldsed,C2) random=REAL(oldsed)/REAL(C2) RETURN END C ************************************************************* C * Uses formulaL: * C * old_cost-cost_of_old_swap+cost_of_new_swap * C * Note: this is fairly generic to objective functions * C * that rely on a linear sum. Can be used in INTSA * C ************************************************************* INTEGER FUNCTION partial_calculate_fitness(fitness,location1 &,location2,city_matrix,solution,old_solution,num_cities) INCLUDE 'SAPARAM.INC' INTEGER solution(max_cols),city_matrix(max_cities &,max_cities),old_solution(max_cols),fitness,location1,location2 &,num_cities,tour_length,previous_city_1,next_city_1, &previous_city_2,next_city_2 IF (location1.EQ.1) THEN previous_city_1=num_cities ELSE previous_city_1=location1-1 ENDIF IF (location1.EQ.num_cities) THEN next_city_1=1 ELSE next_city_1=location1+1 ENDIF IF (location2.EQ.1) THEN previous_city_2=num_cities ELSE previous_city_2=location2-1 ENDIF IF (location2.EQ.num_cities) THEN next_city_2=1 ELSE next_city_2=location2+1 ENDIF c PRINT *,"CC",location1,location2,previous_city_1,previous_city_2 c &,next_city_1,next_city_2,fitness tour_length=fitness-(city_matrix(old_solution(location1), &old_solution(previous_city_1))+city_matrix(old_solution(location2) &,old_solution(next_city_2)))+city_matrix(solution(location1), &solution(previous_city_1))+city_matrix(solution(location2), &solution(next_city_2)) partial_calculate_fitness=tour_length RETURN END C *************************************************************** C * This is a full calculation of the cost function. * C *************************************************************** SUBROUTINE calculate_fitness(solution,city_matrix, &qap_order,fitness) INCLUDE 'SAPARAM.INC' INTEGER solution(max_cols),city_matrix(max_cities &,max_cities) INTEGER gene,sum,sum1,c1,c2,d1,d2,i,j,qap_order,fitness sum1=0 DO i=2,qap_order sum1=sum1+city_matrix(solution(i),solution(i-1)) c PRINT *,sum1,city_matrix(solution(i),solution(i-1)), c & solution(i),solution(i-1),i,i-1 ENDDO sum1=sum1+city_matrix(solution(qap_order),solution(1)) c PRINT *,"sum1=",sum1 fitness=sum1 RETURN END C ************************************************************ C *Prefered random number generator * C ************************************************************ DOUBLE PRECISION FUNCTION RN55() C THIS FUNCTION MUST BE INITIALIZED BY CALLING RNSD(N), C WHERE 0