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. 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

You need to login to post a comment.