Clear Filters
Clear Filters

Solving the hub location problem

46 views (last 30 days)
马清
马清 on 19 Sep 2024 at 15:06
Commented: 马清 on 25 Sep 2024 at 11:11
Problem description:
When solving the Hub Location Problem (HLP), a genetic algorithm was used. Below is the genetic algorithm code for the entire model, but when running the main program `data.m`, the following error message keeps appearing. Despite multiple modifications, the same error persists. Please carefully review all the code for solving the HLP and ensure that the code is corrected to produce the correct results.
Error:
Error using evaluateFitness (line 12)
Individual length does not match the expected total length. Expected: 186, Got: 100
Error in data (line 146)
[obj1, obj2] = evaluateFitness(individual, num_hubs, num_manufacturers, num_customers,
...
data.m (main program file)
beta = 0.7;
num_hubs = 6;
num_manufacturers = 4;
num_customers = 20;
% 1. distance matrix
manufacturers = [113.5494,22.7076; 113.6798,22.6487;
113.3497,22.8446; 113.3735,22.9538]; % Manufacturers' coordinates
hubs = [113.2734,23.1352; 113.3235,23.0990; 113.3685,23.1304;
106.6296,26.6836; 113.2266,23.1875; 113.5936,22.8079]; % Hub centers' coordinates
customers = [113.2702,23.1671; 113.2910,23.2179; 113.2630,23.2053;
113.3101,23.2404; 113.2760,23.0623; 113.2226,23.2340;
113.2947,23.2567; 113.3622,23.2983; 113.4110,23.3917;
113.2777,23.1401; 113.3057,23.1285; 113.2864,23.1507;
113.3144,23.1049; 113.2791,23.0815; 113.2316,23.1262;
113.2508,23.0613; 113.3601,23.1581; 113.3235,23.1685;
113.3934,23.1393; 113.4756,23.1078]; % Customers' coordinates
% Matrix of the distance between the manufacturer and the hub center d_ki
d_ki = [0.5089 0.4519 0.4599 7.9807 0.3643 0.1767;
0.5784 0.1096 0.6339 0.5742 0.2466 7.7066;
0.5735 8.1232 0.7041 0.1810 0.2072 0.2760;
0.3004 0.2557 0.2864 7.7394 0.1536 0.2641];
% The distance matrix between hub centers d_ii
d_ii = [0 0.0618 0.0952 7.5320 0.0702 0.4579;
0.0618 0 0.0549 7.5933 0.1312 0.3971;
0.0952 0.0549 0 7.6183 0.1530 0.3933;
7.5320 7.5933 7.6183 0 7.4661 7.9698;
0.0702 0.1312 0.1530 7.4661 0 0.5280;
0.4579 0.3971 0.3933 7.9698 0.5280 0];
% The distance matrix of customers to the hub center d_ij
d_ij = [0.0321 0.0846 0.0709 0.1114 0.0729 0.1111 0.1234 0.1857 0.2911 0.0065 0.0330 0.0202 0.0510 0.0540 0.0428 0.0773 0.0897 0.0602 0.1201 0.2040;
0.0865 0.1233 0.1223 0.1420 0.0600 0.1685 0.1603 0.2030 0.3055 0.0615 0.0345 0.0636 0.0108 0.0477 0.0958 0.0819 0.0695 0.0695 0.0807 0.1524;
0.1049 0.1169 0.1294 0.1245 0.1149 0.1789 0.1463 0.1680 0.2647 0.0913 0.0628 0.0846 0.0598 0.1019 0.1370 0.1365 0.0289 0.0590 0.0264 0.1095;
7.5142 7.5090 7.4900 7.5156 7.5689 7.4409 7.4945 7.5358 7.5382 7.5335 7.5637 7.5362 7.5825 7.5625 7.4994 7.5473 7.5979 7.5607 7.6362 7.7236;
0.0481 0.0712 0.0405 0.0988 0.1346 0.0467 0.0971 0.1751 0.2751 0.0697 0.0987 0.0702 0.1205 0.1183 0.0615 0.1285 0.1367 0.0987 0.1736 0.2614;
0.4833 0.5096 0.5169 0.5171 0.4069 0.5650 0.5392 0.5423 0.6117 0.4584 0.4309 0.4603 0.4076 0.4169 0.4820 0.4263 0.4209 0.4505 0.3872 0.3223];
% 2. traffic matrix
h_ki = [30 35 60 90 120 100;
25 36 58 70 150 85;
15 20 38 50 95 55;
90 115 160 260 320 190];
g_ii = [0 2 5 7 10 15;
7 0 6 8 12 10;
15 10 0 6 18 12;
6 25 10 0 5 9;
6 8 7 4 0 10;
2 8 10 12 4 0];
s_ij = [3 4 5 8 10 9 2 5 4 6 25 15 2 6 4 5 6 8 4 6;
8 7 12 14 10 12 14 25 18 15 6 7 5 1 5 7 8 9 6 7;
5 6 4 7 6 8 10 14 3 6 2 15 17 15 1 2 5 8 25 10;
3 6 7 9 13 11 1 2 5 7 10 14 5 4 2 6 8 10 3 2;
1 4 5 8 8 7 9 6 10 15 5 4 2 4 6 10 7 4 5 6;
14 12 5 2 3 4 10 9 3 4 3 2 10 15 9 2 3 8 20 12];
% 3. Transportation cost matrix
c_ki = [0.3448 0.3384 0.3393 1.1840 0.3527 0.3000;
0.3589 0.3522 0.3521 1.2000 0.3668 0.3080;
0.3214 0.3164 0.3199 1.1569 0.3286 0.3154;
0.3110 0.3049 0.3075 1.1532 0.3187 0.3174];
c_ii = [0.0000 0.2054 0.2084 0.8615 0.2062 0.2402;
0.2054 0.0000 0.2048 0.8669 0.2115 0.2349;
0.2084 0.2048 0.0000 0.8691 0.2134 0.2345;
0.8615 0.8669 0.8691 0.0000 0.8558 0.9000;
0.2062 0.2115 0.2134 0.8558 0.0000 0.2464;
0.2402 0.2349 0.2345 0.9000 0.2464 0.0000];
c_ij = [0.3030 0.3093 0.3115 1.1756 0.3049 0.3556 0.3091 0.3136 0.3129 1.1750 0.3075 0.3587 0.3075 0.3135 0.3143 1.1728 0.3040 0.3595 0.3122 0.3158;
0.3138 1.1757 0.3108 0.3595 0.3077 0.3062 0.3126 1.1820 0.3149 0.3467 0.3122 0.3189 0.3201 1.1670 0.3047 0.3651 0.3136 0.3179 0.3163 1.1733;
0.3106 0.3621 0.3209 0.3229 0.3188 1.1781 0.3197 0.3625 0.3332 0.3349 0.3301 1.1784 0.3313 0.3706 0.3000 0.3064 0.3099 1.1778 0.3074 0.3527;
0.3031 0.3033 0.3066 1.1814 0.3108 0.3495 0.3016 0.3067 0.3091 1.1781 0.3074 0.3529 0.3052 0.3005 0.3062 1.1835 0.3133 0.3468 0.3055 0.3048;
0.3111 1.1812 0.3130 0.3479 0.3042 0.3104 0.3152 1.1739 0.3064 0.3555 0.3083 0.3088 0.3152 1.1794 0.3142 0.3490 0.3097 0.3073 0.3026 1.1853;
0.3152 0.3483 0.3063 0.3073 0.3061 1.1810 0.3108 0.3495 0.3132 0.3087 0.3023 1.1898 0.3195 0.3444 0.3230 0.3170 0.3120 1.2000 0.3297 0.3368];
% 4. Fixed costs of alternative hub constructionf_i
f_i = [5400 5500 6500 7200 8000 7500];
% 5. Matrix of time
t_ij = [3.2100 5.4650 2.6470 6.8142 4.8100 5.8330 5.4650 3.6330 2.8945 6.6090 4.7120 5.6960 3.8090 5.2230 2.4050 6.4900 4.4050 5.6690 5.8114 3.2020;
2.6450 6.8156 4.9880 5.6710 3.6290 5.4060 2.8149 6.6890 4.6746 5.4650 5.2340 3.8603 2.2463 6.4945 4.8710 5.8923 5.6510 3.4560 2.8445 6.6358;
4.7510 5.4230 3.2911 5.3055 2.6470 6.6382 4.2751 5.6117 3.0650 5.6150 2.9130 6.5335 4.6970 5.4584 3.3300 5.3450 2.6280 6.5637 4.9870 5.4309;
3.2020 5.6360 2.8460 6.5362 4.7020 5.4603 3.5100 5.1080 2.5980 6.5825 4.2050 5.4076 3.5400 5.4770 2.1019 6.5625 4.1830 5.4169 3.4280 5.9580;
2.1370 6.4994 4.6150 5.4820 3.7730 5.8190 2.1365 6.5473 4.2850 5.4263 3.8970 5.6950 2.2890 6.5979 4.3670 5.4209 3.6020 5.6950 2.5900 6.5607;
4.9870 5.4505 3.2010 5.8070 2.2640 6.6362 4.7360 5.3872 3.2040 5.1524 2.1095 6.7236 4.2614 5.3223 3.3015 5.2678 4.3567 6.8574 4.2587 5.2698];
% The maximum tolerable time that the customer can accept T_j
T_j = [4.7965 4.8349 4.6668 4.8555 4.7798 4.7664 4.7936 4.6281 4.7137 4.7141 4.7304 4.5685 4.5469 4.6866 4.6645 4.7111 4.8142 4.6719 4.4622 4.8853];
% 6. Binary variables
% variable x_ki: whether each manufacturer k is connected to each hub i
x_ki = zeros(num_manufacturers, num_hubs); % The matrix initialized to zero
% Variable y i: whether each candidate hub center i is selected
y_i = zeros(num_hubs, 1); % The matrix initialized to zero
% Variable y ii: Whether there is a connection between each pair of hub centers i and i
y_ii = zeros(num_hubs, num_hubs); %The matrix initialized to zero
% Variable z ij: Whether each customer j is served by hub i
z_ij = zeros(num_hubs, num_customers); % The matrix initialized to zero
% parameter setting
popSize = 100; % population size
numGenerations = 150; % iterations
crossoverRate = 0.8; % crossover rate
mutationRate = 0.01; % mutation rate
% Calculate the length of each individual
individual_length = (num_manufacturers * num_hubs) + num_hubs + (num_hubs * num_hubs) + (num_hubs * num_customers);
% Initialize the population (each individual is a binary vector representing the hub center selection)
population = randi([0 1], popSize, individual_length); % Generate a 100x186 0/1 matrix
% The fitness of each generation was recorded
obj1_values = zeros(numGenerations, 1); % Save the value of the optimal obj1 for each generation
obj2_values = zeros(numGenerations, 1); % Save the value of the optimal obj2 for each generation
% 2D floor plan display and fitness change mapping
figure;
hold on;
% Plot the scatter of manufacturers and customers
scatter(manufacturers(:, 1), manufacturers(:, 2), 'b', 'filled'); % manufactures
scatter(customers(:, 1), customers(:, 2), 'g', 'filled'); % customers
xlabel('X axis');
ylabel('Y axis');
title('Distribution of manufacturers and customers');
% Initializes the graphics store
fitnessHistory_1 = zeros(numGenerations, 1); % Used to store the fitness of obj 1 for each generation
fitnessHistory_2 = zeros(numGenerations, 2); % Used to store the fitness of obj 2 for each generation
% Main loop of genetic algorithm
for gen = 1:numGenerations
% For each individual in the population, two objective values are calculated
fitness = zeros(popSize, 2); % Used to store two target values for each individual
for i = 1:popSize
individual = population(i, :); % The ith individual is selected from the population
[obj1, obj2] = evaluateFitness(individual, num_hubs, num_manufacturers, num_customers, ...
d_ki, d_ii, d_ij, h_ki, g_ii, s_ij, ...
c_ki, c_ii, c_ij, f_i, t_ij, T_j);
fitness(i, :) = [obj1, obj2]; % Storage fitness
end
% Save the best fitness value for this generation
fitnessHistory(gen, :) = [min(fitness(:, 1)), max(fitness(:, 2))]; %The fitness values of each generation were recorded
% Select operation (according to fitness)
selected = rouletteWheelSelection(fitness, popSize);
% Crossover operation
offspring = crossover(selected, crossoverRate, num_hubs);
%Mutation operation
offspring = mutation(offspring, mutationRate);
% Renewal of population
population = offspring;
% Draw fitness changes (every two generations to reduce the number of images)
if mod(gen, 2) == 0
subplot(2, 1, 1);
plot(1:gen, fitnessHistory(1:gen, 1), 'b-o', 'DisplayName', 'total cost');
xlabel('iterations');
ylabel('total cost');
title('Change in fitness: total cost');
legend('Location', 'best');
subplot(2, 1, 2);
plot(1:gen, fitnessHistory(1:gen, 2), 'g-o', 'DisplayName', 'customer satisfaction');
xlabel('iterations');
ylabel('customer satisfaction');
title('Change in fitness: customer satisfaction');
legend('Location', 'best');
drawnow; % Dynamic update image
end
end
Error using solution>evaluateFitness (line 194)
Individual length does not match the expected total length. Expected: 186, Got: 100
% Output the selected hub center
selected_hubs = find(population(bestIndexObj1, 1:num_hubs) == 1); % Save the selected hub center number
fprintf('The selected hub center number is:\n');
disp(selected_hubs);
evaluateFitness.m
function [obj1, obj2] = evaluateFitness(individual, num_hubs, num_manufacturers, num_customers, ...
d_ki, d_ii, d_ij, h_ki, g_ii, s_ij, ...
c_ki, c_ii, c_ij, f_i, t_ij, T_j)
% Calculate the expected length in evaluateFitness.m
expected_length = (num_manufacturers * num_hubs) + num_hubs + (num_hubs * num_hubs) + (num_hubs * num_customers);
% Get the length of individual
individual_length = length(individual);
% Length check
if individual_length ~= expected_length
error('Individual length does not match the expected total length. Expected: %d, Got: %d', expected_length, individual_length);
end
% Extract binary variables from individual
x_ki = reshape(individual(1:num_manufacturers * num_hubs), num_manufacturers, num_hubs); % Manufacturer to hub connections
y_i = individual(num_manufacturers * num_hubs + 1:num_manufacturers * num_hubs + num_hubs); % Hub center open status
% Calculate y_ii
start_idx_y_ii = num_manufacturers * num_hubs + num_hubs + 1;
y_ii = reshape(individual(start_idx_y_ii:start_idx_y_ii + (num_hubs * num_hubs) - 1), num_hubs, num_hubs); % Hub-to-hub connections
% Calculate z_ij
start_idx_z_ij = start_idx_y_ii + (num_hubs * num_hubs);
z_ij = reshape(individual(start_idx_z_ij:end), num_hubs, num_customers); % Customer service status
% Initialize objective functions
obj1 = 0; % Objective function 1: Total cost of hubs
obj2 = 0; % Objective function 2: Customer satisfaction
% Initialize penalty
penalty = 0;
% Objective function 1 - Total cost of hubs
for k = 1:num_manufacturers
if sum(x_ki(k, :)) ~= 1 % Check constraint: each manufacturer is connected to exactly one hub
penalty = penalty + 1000; % Add penalty
end
for i = 1:num_hubs
if y_i(i) == 1 && x_ki(k, i) == 1 % Only consider open hubs
obj1 = obj1 + d_ki(k, i) * c_ki(k, i) * h_ki(k, i); % Transportation cost
end
end
end
% Cost of building hub centers
for i = 1:num_hubs
if y_i(i) == 1 % Only calculate fixed cost for open hub centers
obj1 = obj1 + f_i(i);
end
end
% Transportation cost between hub centers
for i = 1:num_hubs
for i_prime = 1:num_hubs
if i ~= i_prime && y_i(i) == 1 && y_ii(i, i_prime) == 1 % Only consider open hubs
obj1 = obj1 + d_ii(i, i_prime) * c_ii(i, i_prime) * g_ii(i, i_prime);
end
end
end
% Transportation cost from hub centers to customers
for i = 1:num_hubs
for j = 1:num_customers
if y_i(i) == 1 && z_ij(i, j) == 1 % Only consider customers served by open hubs
obj1 = obj1 + d_ij(i, j) * c_ij(i, j) * s_ij(i, j);
end
end
end
% Objective function 2 - Customer satisfaction
for i = 1:num_hubs
for j = 1:num_customers
if y_i(i) == 1 && z_ij(i, j) == 1 % Only consider open hubs and their served customers
% Customer satisfaction function
if t_ij(i, j) > 0 && t_ij(i, j) <= 1
satisfaction = 1.00;
elseif t_ij(i, j) > 1 && t_ij(i, j) <= 2
satisfaction = 0.85;
elseif t_ij(i, j) > 2 && t_ij(i, j) <= 3
satisfaction = 0.75;
elseif t_ij(i, j) > 3 && t_ij(i, j) <= 4
satisfaction = 0.45;
elseif t_ij(i, j) > 4 && t_ij(i, j) <= 5
satisfaction = 0.20;
else
satisfaction = 0.00; % Out-of-range time, satisfaction is 0
end
obj2 = obj2 + satisfaction; % Accumulate customer satisfaction
end
end
end
% Check constraints and add penalties
for j = 1:num_customers
if sum(z_ij(:, j)) ~= 1
penalty = penalty + 1000; % Add penalty, no customer should be served by multiple hubs
end
end
for k = 1:num_manufacturers
for i = 1:num_hubs
if x_ki(k, i) > y_i(i)
penalty = penalty + 1000; % Add penalty, manufacturer connected to a non-open hub
end
end
end
for i = 1:num_hubs
for j = 1:num_customers
if z_ij(i, j) > y_i(i)
penalty = penalty + 1000; % Add penalty, customer should not be served by a non-open hub
end
end
end
% Flow balance constraint
for i = 1:num_hubs
if sum(h_ki(:, i) .* x_ki(:, i) .* y_i(i)) ~= sum(s_ij(i, :) .* z_ij(i, :) .* y_i(i))
penalty = penalty + 1000; % Add penalty, flow balance mismatch
end
end
% Maximum allowable time for customers
for i = 1:num_hubs
for j = 1:num_customers
if t_ij(i, j) * z_ij(i, j) * y_i(i) > T_j(j)
penalty = penalty + 1000; % Add penalty, customer service time exceeds maximum allowable time
end
end
end
% Symmetry constraint
for i = 1:num_hubs
for i_prime = 1:num_hubs
if y_ii(i, i_prime) ~= y_ii(i_prime, i)
penalty = penalty + 1000; % Add penalty, symmetry mismatch
end
end
end
% Constraint 7: y_ii <= min(y_i, y_i')
for i = 1:num_hubs
for i_prime = 1:num_hubs
if y_ii(i, i_prime) > min(y_i(i), y_i(i_prime))
penalty = penalty + 1000; % Add penalty, constraint mismatch
end
end
end
% Ensure all binary variables are 0 or 1
if any([x_ki(:); y_i(:); z_ij(:); y_ii(:)] ~= round([x_ki(:); y_i(:); z_ij(:); y_ii(:)]))
penalty = penalty + 1000; % Add penalty
end
% Add penalty to obj1 and obj2
obj1 = obj1 + penalty;
obj2 = obj2 - penalty; % Adjust the impact of penalty on obj2 as needed
end
generateIndividual.m
function individual = generateIndividual(num_manufacturers, num_hubs, num_customers)
% Calculate the expected length of the individual
expected_length = (num_manufacturers * num_hubs) + num_hubs + (num_hubs * num_hubs) + (num_hubs * num_customers);
% Generate the individual
individual = zeros(expected_length, 1);
% Fill in each part
individual(1:num_manufacturers * num_hubs) = randi([0, 1], num_manufacturers * num_hubs, 1); % Manufacturer to hub connections
individual(num_manufacturers * num_hubs + 1:num_manufacturers * num_hubs + num_hubs) = randi([0, 1], num_hubs, 1); % Hub selection
individual(num_manufacturers * num_hubs + num_hubs + 1:num_manufacturers * num_hubs + num_hubs + num_hubs * num_hubs) = randi([0, 1], num_hubs * num_hubs, 1); % Connections between hubs
individual(num_manufacturers * num_hubs + num_hubs + num_hubs * num_hubs + 1:end) = randi([0, 1], num_hubs * num_customers, 1); % Connections from hubs to customers
end
selection.m
% Selection function
function selected = selection(pop, obj1, obj2, popSize, alpha)
% alpha is a parameter controlling the weight between objective function 1 and objective function 2, where 0 <= alpha <= 1
% Calculate the combined fitness of each individual
fitness = alpha * obj1 + (1 - alpha) * obj2;
% Roulette wheel selection
prob = exp(fitness - max(fitness)); % Prevent numerical overflow
prob = prob / sum(prob);
selected = pop(randsample(size(pop, 1), popSize, true, prob), :);
end
rouletteWheelSelection.m
function selected = rouletteWheelSelection(fitness, popSize)
% Calculate the total fitness of all individuals
totalFitness = sum(fitness);
% Calculate the cumulative probability distribution
cumProb = cumsum(fitness) / totalFitness;
% Initialize the selection result
selected = zeros(1, popSize);
for i = 1:popSize
r = rand(); % Generate a random number in the interval [0, 1]
index = find(cumProb >= r, 1); % Find the first individual with cumulative probability >= r
if ~isempty(index)
selected(i) = index;
else
% Handle the case where no suitable individual is found
selected(i) = length(cumProb); % Choose the last individual as a default value
end
end
end
crossover.m
function offspring = crossover(selected, crossoverRate, num_hubs)
numSelected = size(selected, 1); % Get the number of selected individuals
numGenes = size(selected, 2); % Get the number of genes
popSize = size(selected, 1); % Get the population size
offspring = zeros(popSize, numGenes); % Initialize the offspring matrix with size popSize x numGenes
% Perform crossover based on the crossover rate
for i = 1:2:popSize
if rand <= crossoverRate
% Randomly select parents
parent1 = selected(randi(numSelected), :);
parent2 = selected(randi(numSelected), :);
% Randomly choose a crossover point
crossoverPoint = randi(numGenes - 1);
% Generate offspring
offspring(i, :) = [parent1(1:crossoverPoint), parent2(crossoverPoint+1:end)];
if i < popSize
offspring(i+1, :) = [parent2(1:crossoverPoint), parent1(crossoverPoint+1:end)];
end
else
% If no crossover, copy directly
offspring(i, :) = selected(i, :);
if i < popSize
offspring(i+1, :) = selected(i+1, :);
end
end
end
end
mutation.m
function mutated = mutation(offspring, mutationRate)
numOffspring = size(offspring, 1);
numGenes = size(offspring, 2);
mutated = offspring;
for i = 1:numOffspring
for j = 1:numGenes
if rand < mutationRate
mutated(i, j) = 1 - mutated(i, j); % Flip the bit
end
end
end
end

Accepted Answer

Torsten
Torsten on 19 Sep 2024 at 19:12
Moved: Torsten on 19 Sep 2024 at 19:12
You set
% Renewal of population
population = offspring;
and offspring has only size 1x100 instead of 1x186.
Thus in your loop
% Main loop of genetic algorithm
for gen = 1:numGenerations
% For each individual in the population, two objective values are calculated
fitness = zeros(popSize, 2); % Used to store two target values for each individual
for i = 1:popSize
individual = population(i, :); % The ith individual is selected from the population
[obj1, obj2] = evaluateFitness(individual, num_hubs, num_manufacturers, num_customers, ...
d_ki, d_ii, d_ij, h_ki, g_ii, s_ij, ...
c_ki, c_ii, c_ij, f_i, t_ij, T_j);
fitness(i, :) = [obj1, obj2]; % Storage fitness
end
when gen = 2, population(i, :) has changed size, and evaluateFitness errors.
  3 Comments
Torsten
Torsten on 20 Sep 2024 at 8:06
I have no knowledge about the background of the hub location problem.
If you give a mathematical description of what you are trying to solve, maybe someone is able to help. Else I would contact the authors of the code and ask what went wrong.
马清
马清 on 25 Sep 2024 at 11:11
I'm so sorry for responding to your valuable answer about solving the HLP so late, as I have had to handle many issues regarding my company in recent days. Since the previous HLP file is not related to my research field, I need to modify the HLP file. In the new file, I have added specific details about the solved HLP (e.g., the background of the problem, some relevant assumptions) so that readers can better understand and solve the HLP. Additionally, I still hope to use the genetic algorithm (GA) to solve the model. Therefore, could you please provide me with detailed MATLAB programming for GA to handle the HLP model? I would be extremely grateful! I'm looking forward to your response.

Sign in to comment.

More Answers (0)

Categories

Find more on Sensitivity Analysis in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!