function [routes, totalCost] = cvrp_saving(d, q)
% cvrp_saving 使用节约算法求解cvrp问题
% d: 距离矩阵，描述任意两个客户之间的距离，其中d(i,j)表示客户i和客户j之间的距离。
% q: 向量，表示每个客户的需求量。
% routes: 返回所有车辆的行驶路径，其中routes{k}表示k号车的行驶路径。
% totalCost: 所有车总共耗时。

% 初始化变量
n = length(q); % 客户数目 
capacity = sum(q); % 所有客户总需求量 

% 节省域初始化：将所有不合理的对来保存在一个二位数组中 
S = zeros(n, n);

% 构建Saving Domain S (n x n) 
for i=1:(n-1)      % 第一个for循环用于遍历所有不合理对     

    for j=(i+1):n      % 第二个for循环用于遍历当前不合理对中剩余部分     

        if q(i)+q(j) <= capacity && d(i,j)-d(0,i)-d(0,j) > 0       % 当前不合理对是否可以加入Saving Domain         

            S(i, j)= d(i, j)-d (0 , i ) - d (0 , j );           % 计算当前不合理对的saving                                   value    

        end    

    end    

end    

while ~isempty([find (S==max ( max (S)))])          % 循坏：当Saving Domain仍然存在时              [m , n] = find (S==max ( max (S)));           sMax=max ( max (S));               S[m , n]= 0 ;           q[m]=q[m]+q[n];            q[n]= 0 ;       end         routes={};         routes{1}=[0 m];         k= 1 ;             while ~isempty ([find ((q!= 0 ) & ((sum((~ismember ([0 m], routes{k})))!= 0 ))])              [a , b] = find ((q!= 0 ) & ((sum((~ismember ([0 m], routes{k})))!= 0 )));              minDij=- 1 ;              minJ=- 1 ;                  for j=(b+ 1 ):length([find ((q!= 0 ) & ((sum((~ismember ([0 m], routes{k})))!= 0 )))])                  if minDij == - 1 || d([a b]) < minDij                      minDij = d([a b]);                      minJ=[b a];                  end             end             if ~isequal ([minJ],[ - 1 -1 ])                 k++;                 routes{k}=[minJ];             else                 break;             end         end          totalCost= sum((cellfun (@length ,routes))*capacity)* 10 / 3600;