Posted By

mertnuhoglu on 12/07/07


Tagged

matlab heuristics simulated-annealing location-allocation


Versions (?)

Who likes this?

1 person have marked this snippet as a favorite

mertnuhoglu


Simulated annealing for location allocation problem


 / Published in: MatLab
 

Simulated annealing heuristics applied to solve continuous location allocation problem.

  1. function [report,result] = solve_problems()
  2. data_files={'bon287';'p654'};
  3. p=[2 4 5 6];
  4. n_runs=5;
  5. k=0;
  6. for run=1:n_runs
  7. for i=1:length(data_files)
  8. for j=1:length(p)
  9. k=k+1;
  10. data_file=char(data_files(i));
  11. locations=load([data_file '_locations.txt']);
  12. demands=load([data_file '_demands.txt']);
  13. [z,x,cycles,debug]=solve_location_allocation(locations,demands,p(j));
  14. dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_cycles.txt'],cycles);
  15. dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_debug.txt'],debug);
  16. dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_result.txt'],[z x]);
  17. result(run,i,j)=z;
  18. end
  19. end
  20. end
  21.  
  22. function [z_best,x_best,cycles,debug]=solve_location_allocation(custLocs,demands,p,acceptanceRule)
  23. if nargin < 4
  24. end
  25. n_c=length(demands); % number of customers
  26. tc=0; % termination counter
  27. ip=0.3; % updates per cycle bound for temperature
  28. r=0.9; % cooling ratio
  29. fp=0.02; % updates per cycle bound for termination
  30. gamma=1.1; % cycle growth factor
  31. p_1=0.95; % initial probability
  32. k=1; % number of changes in the neighborhood function
  33. thresholdParameter=0.1;
  34. T=initializeTemperature(p_1,p,k,custLocs,demands);
  35. ns=nchoosek(n_c,k)*(p-1)^k;
  36. L=2*ns;
  37. x=generateArbitrarySolution(n_c,p); % customer/facility assignments
  38. z_best=inf;
  39. cycle=1;
  40. cycles(1,:)=[z_best L T tc]; % parameters stored for debugging per cycle
  41. debug(1,:)=[z_best inf inf 0 L 0]; % parameters stored for debugging (per iteration)
  42. while(~stoppingCondition(tc,cycles))
  43. j=0;
  44. for i=1:L
  45. facLocs=single_facility_optimization(x,custLocs,demands);
  46. x=findAssignments(facLocs,custLocs);
  47. z=f(x,demands,facLocs,custLocs);
  48. x_=pickANeighbor(x,k,p);
  49. facLocs_=single_facility_optimization(x_,custLocs,demands);
  50. x_=findAssignments(facLocs_,custLocs);
  51. z_=f(x_,demands,facLocs_,custLocs);
  52. delta=z-z_;
  53. if(delta<0)
  54. x=x_;
  55. j=j+1;
  56. if(z_<z_best)
  57. i
  58. z_best=z_
  59. x_best=x_;
  60. facLocs_best=facLocs_;
  61. end
  62. else
  63. if(acceptanceRule(delta,T,z_,z,thresholdParameter))
  64. x=x_;
  65. j=j+1;
  66. end
  67. end
  68. debug(end+1,:)=[z_best z z_ j L i];
  69. end
  70. tc=tc+changeTerminationCounter(j,L,fp);
  71. T=decreaseTemperature(j,L,ip,r,T);
  72. thresholdParameter=thresholdParameter*r;
  73. L=round(L*gamma)
  74. cycle=cycle+1
  75. cycles(cycle,:)=[z_best L T tc];
  76. end
  77.  
  78. function result=thresholdAcceptance(delta,T,z_,z,mu)
  79. result=z_<=(1+mu)*z;
  80.  
  81. function result=probabilisticAcceptance(delta,T,z_,z,mu)
  82. result=exp(-delta/T)>rand;
  83.  
  84. function result=stoppingCondition(terminationCounter,cycles)
  85. % result = terminationCounter >= 5; % alternative stopping condition
  86. z_best=cycles(:,1);
  87. if(length(z_best)<2)
  88. result=false;
  89. return;
  90. end
  91. change=(z_best(end)-z_best(end-1))/z_best(end);
  92. eps_1=0.03;
  93. if(change < eps_1)
  94. result=true;
  95. else
  96. result=false;
  97. end
  98.  
  99. function result=changeTerminationCounter(j,L,fp)
  100. if(j/L <= fp)
  101. result=1;
  102. else
  103. result=0;
  104. end
  105.  
  106. function T=decreaseTemperature(j,L,ip,r,T)
  107. if(j/L > ip)
  108. T=T/2;
  109. else
  110. T=r*T;
  111. end
  112.  
  113. function result=distance(x,y,degree)
  114. if nargin < 3
  115. degree=2; % default: euclidean distance
  116. end
  117. result=(abs(x(1)-y(1))^degree+abs(x(2)-y(2))^degree)^(1/degree);
  118.  
  119. % neighborhood structure is based on the customer-facility assignments
  120. % pick a neighbor of the current solution x where k is the number of changing assignments, p is the number of facilities
  121. function x_=pickANeighbor(x,k,p)
  122. changingCustomers=unidrnd(length(x),1,k);
  123. x_=x;
  124. for i=1:k
  125. newFacility=unidrnd(p-1);
  126. oldFacility=x(changingCustomers(i));
  127. if (newFacility >= oldFacility)
  128. newFacility=newFacility+1;
  129. end
  130. x_(changingCustomers(i))=newFacility;
  131. end
  132.  
  133. function T=initializeTemperature(p_1,p,k,custLocs,demands)
  134. n=100;
  135. n_c=length(demands);
  136. for i=1:n
  137. x=generateArbitrarySolution(n_c,p); % customer-facility assignments
  138. x_=pickANeighbor(x,k,p);
  139. facLocs=single_facility_optimization(x,custLocs,demands);
  140. facLocs_=single_facility_optimization(x_,custLocs,demands);
  141. Delta(i)=abs(f(x,demands,facLocs,custLocs)-f(x_,demands,facLocs_,custLocs));
  142. end
  143. T=mean(Delta)/log(1/p_1);
  144.  
  145. function x=generateArbitrarySolution(n_c,p)
  146. x=unidrnd(p,1,n_c);
  147.  
  148. function result=findAssignments(facLocs,custLocs)
  149. n_c=length(custLocs);
  150. p=length(facLocs);
  151. distances=distancesFromFacilities(facLocs,custLocs);
  152. [minD,result]=min(distances');
  153.  
  154. % objective function value of a solution x
  155. function result=f(x,demands,facLocs,custLocs)
  156. n_c=length(custLocs);
  157. result=0;
  158. for cust=1:n_c
  159. facility=x(cust);
  160. dist=distance(facLocs(facility,:),custLocs(cust,:));
  161. result=result+dist*demands(cust);
  162. end
  163.  
  164. function result=distancesFromSingleFacility(facilityLocation,customerLocations,degree)
  165. if nargin < 3
  166. degree=2; % default: euclidean distance
  167. end
  168. for cust=1:length(customerLocations)
  169. y=customerLocations(cust,1:2);
  170. result(cust)=distance(facilityLocation,y,degree);
  171. end
  172.  
  173. function distances=distancesFromFacilities(facLocs,custLocs)
  174. n_c=length(custLocs);
  175. for fac=1:length(facLocs)
  176. distances(1:n_c,fac)=distancesFromSingleFacility(facLocs(fac,:),custLocs);
  177. end

Report this snippet  

Comments

RSS Icon Subscribe to comments
Posted By: wst42 on January 8, 2008

Well, it would be great if you can post your data file or tell their format

Thanks

Posted By: neb_sar on August 21, 2008

kardeş paylaşım için çok sağol yanlız swt42 nin dediği gibi data dosyalarının ne olduklarını anlatan bir not yazarsan ya da data dosyalarını da bizimle paylaşırsan ya da e mailime gönderirsen çok sevinirim. location problemi ile ilgili tez yapacağım da bu kodlar gerçekten işime yarayabilir. teşekkürler.

Posted By: mertnuhoglu on September 30, 2008

Sorry for late response. snipplr doesn't send email by default when new comments are added. Therefore I wasn't aware of these requests. I put all the data files under http://sites.google.com/site/mertnuhoglu/programming/files/ie517-heuristics

Posted By: anandsinghal on March 27, 2009

singlefacilityoptimization is not defined in the above file,please help me to sort it out as quick as possible ,my project is in between

You need to login to post a comment.