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