请用禁忌搜索算法解决算例，用matlab代码

要使用禁忌搜索算法（Tabu Search, TS）解决给定的物流系统车辆路径优化问题，并提供 MATLAB 代码实现，以下是详细的步骤和代码：

### 1. 问题描述
- **分拨中心**：P (35, 35)
- **客户点**：30个客户点，每个客户点的位置、需求量和服务时间窗口。
- **车辆参数**：
- 车辆数：20辆
- 额定载重量：2000 kg
- 有效服务时长：3.5小时
- 折损成本：50元/小时
- 早到惩罚成本：30元/小时
- 晚到惩罚成本：120元/小时
- 服务时长：10分钟

### 2. 禁忌搜索算法步骤
1. **初始化**：
- 生成初始解（每辆车只服务一个客户点）
- 设置禁忌表长度
- 设置最大迭代次数
2. **邻域搜索**：
- 生成当前解的邻域解
- 选择最佳邻域解
- 更新禁忌表
3. **终止条件**：
- 达到最大迭代次数或最优解长时间未更新

### 3. MATLAB 代码实现

function [bestSolution, bestCost] = tabuSearchVRP()
    % 参数设置
    maxIter = 2000;
    tabuLength = 20;
    numVehicles = 20;
    depot = [35, 35];
    customers = [
        41, 49, 180, 15, 25;
        55, 45, 200, 20, 30;
        10, 43, 300, 35, 45;
        55, 60, 450, 30, 40;
        30, 60, 200, 150, 160;
        50, 35, 100, 150, 160;
        15, 10, 500, 110, 120;
        30, 5, 200, 135, 145;
        5, 80, 200, 50, 60;
        45, 65, 400, 40, 50;
        13, 65, 350, 80, 90;
        2, 60, 600, 100, 110;
        75, 5, 250, 100, 110;
        96, 38, 150, 70, 80;
        36, 26, 150, 170, 180;
        85, 77, 200, 65, 75;
        87, 69, 260, 80, 90;
        53, 3, 500, 130, 140;
        81, 52, 380, 105, 115;
        16, 18, 280, 90, 100;
        4, 18, 100, 70, 80;
        76, 92, 300, 90, 100;
        80, 26, 450, 100, 110;
        93, 8, 300, 70, 80;
        29, 97, 300, 105, 115;
        53, 87, 200, 72, 82;
        67, 34, 200, 125, 135;
        71, 78, 250, 65, 75;
        89, 87, 350, 110, 120;
        53, 30, 550, 13, 23
    ];

    % 计算距离矩阵
    distMatrix = calculateDistanceMatrix(depot, customers);

    % 初始化解
    initialSolution = initializeSolution(numVehicles, size(customers, 1));
    currentSolution = initialSolution;
    bestSolution = initialSolution;
    bestCost = calculateTotalCost(currentSolution, distMatrix, customers);

    % 禁忌表
    tabuList = cell(tabuLength, 1);
    for iter = 1:maxIter
        % 生成邻域解
        neighborhood = generateNeighborhood(currentSolution);
        
        % 选择最佳邻域解
        [newSolution, newCost] = selectBestNeighbor(neighborhood, distMatrix, customers, tabuList);
        
        % 更新当前解
        if newCost < bestCost
            bestSolution = newSolution;
            bestCost = newCost;
        end
        
        currentSolution = newSolution;
        
        % 更新禁忌表
        updateTabuList(tabuList, newSolution);
        
        % 终止条件
        if mod(iter, 100) == 0
            fprintf('Iteration %d: Best Cost = %.2f\n', iter, bestCost);
        end
    end
end

function distMatrix = calculateDistanceMatrix(depot, customers)
    n = size(customers, 1);
    distMatrix = zeros(n+1, n+1);
    
    % 计算分拨中心到所有客户点的距离
    for i = 1:n
        distMatrix(1, i+1) = norm(depot - customers(i, :));
        distMatrix(i+1, 1) = distMatrix(1, i+1);
    end
    
    % 计算客户点之间的距离
    for i = 1:n
        for j = 1:n
            if i ~= j
                distMatrix(i+1, j+1) = norm(customers(i, :) - customers(j, :));
            end
        end
    end
end

function solution = initializeSolution(numVehicles, numCustomers)
    solution = cell(numVehicles, 1);
    for i = 1:numVehicles
        if i <= numCustomers
            solution{i} = i;
        else
            solution{i} = [];
        end
    end
end

function neighborhood = generateNeighborhood(solution)
    numVehicles = length(solution);
    neighborhood = {};
    
    for i = 1:numVehicles
        for j = i+1:numVehicles
            if ~isempty(solution{i}) && ~isempty(solution{j})
                % 交换客户点
                temp = solution{i};
                solution{i} = solution{j};
                solution{j} = temp;
                neighborhood{end+1} = solution;
                
                % 插入客户点
                for k = 1:length(solution{i})
                    temp = solution{i};
                    solution{i}(k) = solution{j}(1);
                    solution{j}(1) = [];
                    neighborhood{end+1} = solution;
                    solution{i} = temp;
                end
                
                for k = 1:length(solution{j})
                    temp = solution{j};
                    solution{j}(k) = solution{i}(1);
                    solution{i}(1) = [];
                    neighborhood{end+1} = solution;
                    solution{j} = temp;
                end
            end
        end
    end
end

function [bestNeighbor, bestCost] = selectBestNeighbor(neighborhood, distMatrix, customers, tabuList)
    bestCost = inf;
    bestNeighbor = [];
    
    for i = 1:length(neighborhood)
        cost = calculateTotalCost(neighborhood{i}, distMatrix, customers);
        if cost < bestCost && ~isInTabuList(neighborhood{i}, tabuList)
            bestCost = cost;
            bestNeighbor = neighborhood{i};
        end
    end
end

function isInTabu = isInTabuList(solution, tabuList)
    isInTabu = false;
    for i = 1:length(tabuList)
        if isequal(solution, tabuList{i})
            isInTabu = true;
            return;
        end
    end
end

function updateTabuList(tabuList, newSolution)
    tabuList = circshift(tabuList, 1);
    tabuList{1} = newSolution;
end

function totalCost = calculateTotalCost(solution, distMatrix, customers)
    totalCost = 0;
    numVehicles = length(solution);
    
    for i = 1:numVehicles
        if ~isempty(solution{i})
            route = [1, solution{i} + 1]; % 加上分拨中心
            routeCost = 0;
            
            for j = 1:length(route)-1
                routeCost = routeCost + distMatrix(route(j), route(j+1));
            end
            
            routeCost = routeCost * 5; % 运输成本
            routeCost = routeCost + 80; % 启动成本
            
            % 计算时间和惩罚成本
            currentTime = 0;
            for j = 2:length(route)-1
                currentTime = currentTime + distMatrix(route(j-1), route(j)) / 60; % 行驶时间
                currentTime = currentTime + 10 / 60; % 服务时间
                startTime = customers(route(j)-1, 3) / 60;
                endTime = customers(route(j)-1, 4) / 60;
                
                if currentTime < startTime
                    routeCost = routeCost + (startTime - currentTime) * 30; % 早到惩罚
                elseif currentTime > endTime
                    routeCost = routeCost + (currentTime - endTime) * 120; % 晚到惩罚
                end
            end
            
            currentTime = currentTime + distMatrix(route(end-1), route(end)) / 60; % 返回分拨中心时间
            if currentTime > 3.5
                routeCost = routeCost + (currentTime - 3.5) * 50; % 超时折损成本
            end
            
            totalCost = totalCost + routeCost;
        end
    end
end

### 4. 运行代码
将上述代码保存为 `tabuSearchVRP.m` 文件，然后在 MATLAB 中运行 `tabuSearchVRP` 函数即可得到最优解和最小成本。

### 5. 结果分析
- **最优解**：输出的最佳配送方案和总成本。
- **性能评估**：可以通过调整参数（如最大迭代次数、禁忌表长度等）来进一步优化算法性能。

希望这段代码能帮助你解决物流系统车辆路径优化问题！如果有任何问题或需要进一步的帮助，请随时联系我。