function [report,result] = solve_problems()
data_files={'bon287';'p654'};
p=[2 4 5 6];
n_runs=5;
k=0;
for run=1:n_runs
for i=1:length(data_files)
for j=1:length(p)
k=k+1;
data_file=char(data_files(i));
locations=load([data_file '_locations.txt']);
demands=load([data_file '_demands.txt']);
[z,x,cycles,debug]=solve_location_allocation(locations,demands,p(j));
dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_cycles.txt'],cycles);
dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_debug.txt'],debug);
dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_result.txt'],[z x]);
result(run,i,j)=z;
end
end
end
function [z_best,x_best,cycles,debug]=solve_location_allocation(custLocs,demands,p,acceptanceRule)
if nargin < 4
acceptanceRule=@probabilisticAcceptance;
end
n_c=length(demands); % number of customers
tc=0; % termination counter
ip=0.3; % updates per cycle bound for temperature
r=0.9; % cooling ratio
fp=0.02; % updates per cycle bound for termination
gamma=1.1; % cycle growth factor
p_1=0.95; % initial probability
k=1; % number of changes in the neighborhood function
thresholdParameter=0.1;
T=initializeTemperature(p_1,p,k,custLocs,demands);
ns=nchoosek(n_c,k)*(p-1)^k;
L=2*ns;
x=generateArbitrarySolution(n_c,p); % customer/facility assignments
z_best=inf;
cycle=1;
cycles(1,:)=[z_best L T tc]; % parameters stored for debugging per cycle
debug(1,:)=[z_best inf inf 0 L 0]; % parameters stored for debugging (per iteration)
while(~stoppingCondition(tc,cycles))
j=0;
for i=1:L
facLocs=single_facility_optimization(x,custLocs,demands);
x=findAssignments(facLocs,custLocs);
z=f(x,demands,facLocs,custLocs);
x_=pickANeighbor(x,k,p);
facLocs_=single_facility_optimization(x_,custLocs,demands);
x_=findAssignments(facLocs_,custLocs);
z_=f(x_,demands,facLocs_,custLocs);
delta=z-z_;
if(delta<0)
x=x_;
j=j+1;
if(z_<z_best)
i
z_best=z_
x_best=x_;
facLocs_best=facLocs_;
end
else
if(acceptanceRule(delta,T,z_,z,thresholdParameter))
x=x_;
j=j+1;
end
end
debug(end+1,:)=[z_best z z_ j L i];
end
tc=tc+changeTerminationCounter(j,L,fp);
T=decreaseTemperature(j,L,ip,r,T);
thresholdParameter=thresholdParameter*r;
L=round(L*gamma)
cycle=cycle+1
cycles(cycle,:)=[z_best L T tc];
end
function result=thresholdAcceptance(delta,T,z_,z,mu)
result=z_<=(1+mu)*z;
function result=probabilisticAcceptance(delta,T,z_,z,mu)
result=exp(-delta/T)>rand;
function result=stoppingCondition(terminationCounter,cycles)
% result = terminationCounter >= 5; % alternative stopping condition
z_best=cycles(:,1);
if(length(z_best)<2)
result=false;
return;
end
change=(z_best(end)-z_best(end-1))/z_best(end);
eps_1=0.03;
if(change < eps_1)
result=true;
else
result=false;
end
function result=changeTerminationCounter(j,L,fp)
if(j/L <= fp)
result=1;
else
result=0;
end
function T=decreaseTemperature(j,L,ip,r,T)
if(j/L > ip)
T=T/2;
else
T=r*T;
end
function result=distance(x,y,degree)
if nargin < 3
degree=2; % default: euclidean distance
end
result=(abs(x(1)-y(1))^degree+abs(x(2)-y(2))^degree)^(1/degree);
% neighborhood structure is based on the customer-facility assignments
% pick a neighbor of the current solution x where k is the number of changing assignments, p is the number of facilities
function x_=pickANeighbor(x,k,p)
changingCustomers=unidrnd(length(x),1,k);
x_=x;
for i=1:k
newFacility=unidrnd(p-1);
oldFacility=x(changingCustomers(i));
if (newFacility >= oldFacility)
newFacility=newFacility+1;
end
x_(changingCustomers(i))=newFacility;
end
function T=initializeTemperature(p_1,p,k,custLocs,demands)
n=100;
n_c=length(demands);
for i=1:n
x=generateArbitrarySolution(n_c,p); % customer-facility assignments
x_=pickANeighbor(x,k,p);
facLocs=single_facility_optimization(x,custLocs,demands);
facLocs_=single_facility_optimization(x_,custLocs,demands);
Delta(i)=abs(f(x,demands,facLocs,custLocs)-f(x_,demands,facLocs_,custLocs));
end
T=mean(Delta)/log(1/p_1);
function x=generateArbitrarySolution(n_c,p)
x=unidrnd(p,1,n_c);
function result=findAssignments(facLocs,custLocs)
n_c=length(custLocs);
p=length(facLocs);
distances=distancesFromFacilities(facLocs,custLocs);
[minD,result]=min(distances');
% objective function value of a solution x
function result=f(x,demands,facLocs,custLocs)
n_c=length(custLocs);
result=0;
for cust=1:n_c
facility=x(cust);
dist=distance(facLocs(facility,:),custLocs(cust,:));
result=result+dist*demands(cust);
end
function result=distancesFromSingleFacility(facilityLocation,customerLocations,degree)
if nargin < 3
degree=2; % default: euclidean distance
end
for cust=1:length(customerLocations)
y=customerLocations(cust,1:2);
result(cust)=distance(facilityLocation,y,degree);
end
function distances=distancesFromFacilities(facLocs,custLocs)
n_c=length(custLocs);
for fac=1:length(facLocs)
distances(1:n_c,fac)=distancesFromSingleFacility(facLocs(fac,:),custLocs);
end