Simulated annealing for location allocation problem


/ Published in: MatLab
Save to your folder(s)

Simulated annealing heuristics applied to solve continuous location allocation problem.


Copy this code and paste it in your HTML
  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

You need to login to post a comment.