www.impianti.dii.unipg.it

Matlab Code for the Saving Matrix Method 

 

 

Note

You can calculate the savings matrix from the distance matrix in one line of Matlab Code

if you name distance the matrix of distances (n,n) ,where all the distances between n-1 points and DC are shown,

the savings matrix is:

 

savings=tril(distance(2:n,1)*ones(1,n-1)+distance(2:n,1)*ones(1,n-1)'-distance(2:n,2:n),-1)

 

In the following find the matlab code for implementing the Matrix Savings Method

The method work with any number of points. You must change the input files.

You can find an example of input files here.

This code integrated this , by adding the capacitated contstraints.

 

function [giri, mezzi] = tsp_cw()
clear 
clc
close all

%apro e chiudo il file dove scrivere i dati in modo da cancellarne il
%precedente contenuto prima di scrivere i dati della nuova simulazione
route=fopen('route10.txt', 'wt');
fclose(route);

%leggo il numero di nodi totali, comprensivi anche del distribution center (DC).
%potrei anche inserirlo da schermo!!! 
%n= input('inserisci il numero di nodi (compreso il distribution center)  ')

num=fopen ('num10.txt', 'r');
n = fscanf(num, '%g ', 1);
numeri=[1:n];          %vettore necessario per l'ordinamento degli elementi
%leggo da file la matrice delle distanze
%_____1 2 3 4 ...
%___1 
%___2 
%___3 
%___4 
%___... 
%-------------ESEMPIO----------
%0 10 17 25 20 15 16
%10 0 11 15 12 7 16
%17 11 0 9 8 12 24
%25 15 9 0 6 11 24
%20 12 8 6 0 7 18
%15 7 12 11 7 0 13
%16 16 24 24 18 13 0
%-------------------------------

dist=fopen ('data10.txt', 'r');
data = fscanf(dist, '%g \n', [n n]);

%leggo da file il quantitativo da distribuire ad ogni nodo (il primo è
%relativo al distribution center --> quantitativo nullo)
cap=fopen ('capacity10.txt', 'r');
capacity = fscanf(cap, '%g \n', [1 n]);

%leggo da file le capacità totali dei mezzi utilizzati
tr=fopen ('truck10.txt', 'r');
trucks = fscanf(tr, '%g \n', [1 n]);

n_trucks= size(trucks,2);  %numero di righe di trucks -> numero mezzi
n_point= n-1;              %numero dei nodi escluso il DC

while size(data,1)>1 &&  n_trucks>0 && n_point>1  
    
[tour,lenght,out]= clarke_wright(data, capacity, trucks(1),n,numeri) %chiamo la function per costruire il giro

%a=size(tour,2);
%f=size(out,1) ;

    if size(out,1)>2  %serve per verificare quanti nodi rimangono non raggiunti e costruire il nuovo giro di consegna per tali punti
   
        if size(tour,2)>=2  %verifico se nel vettore tour ho 2 o + nodi.
           tour = sort(tour, 'descend');  %ordino i nodi 
            
%data(tour(2):tour(end), :)=[]
           for i=1:size(tour,2) -1  %costruisco la nuova matrice data togliendo le righe dei punti già raggiunti dal DC
               data(tour(i),:)=[]
           end
    
%data(:,tour(2):tour(end))=[]
           for j=1:size(tour,2)-1   %costruisco la nuova matrice data togliendo le colonne dei punti già raggiunti dal DC
               data(:,tour(j))=[]
               capacity(:,tour(j))=[]  %costruisco il nuovo vettore capacity per aggiornare le quantità
               numeri(tour(j))=[]
           end
      
        end
   
n_point= size(out,1)  %aggiorno il numero di nodi
%size_tr=size(trucks,2);

        if size(trucks,2)>1    %aggiorno il numero di mezzi 
           trucks(:,1)=[]
           n_trucks=n_trucks-1  
        else 
           n_trucks=0
        end

n=size(out,1)+1  % i nodi che ho sono quelli rimasti non serviti nel giro precedente + il DC.

    else   %se ho 2 nodi non raggiunti dal giro precedente.verifico le capacità.
        
        if size(out,1)==2 && n_trucks>1   %ho 2 nodi non raggiunti nel giro precedente 
           route=fopen('route10.txt', 'at');
           girofattibile = 0;
           ifatto=0;
           i1fatto=0;
           
         for i=2:n_trucks    
             
            if (capacity(out(1)) + capacity (out(2))) <= trucks(i) %somma delle 2 capacità < capacità del mezzo -> il mezzo può servire i 2 nodi.
                fprintf(route, '\n %s \n', 'Il giro è formato dai punti out(1) e out(2)');
                girofattibile = 1;
            break
            end  
          %somma delle 2 capacità > capacità del mezzo -> il mezzo NON può servire i 2 nodi.
          fprintf(route, '\n %s', 'Non si può fare il giro con il truck! La somma dei quantitativi dei 2 nodi eccede la disponibilità del truck che ha capacità:  '); 
          fprintf(route, '%d \n', trucks(i)); 
       
         end
         
         if girofattibile == 0
            
                 if capacity(out(1))< capacity(out(2))
                    
                 for i=2:n_trucks
                 
                     if trucks(i) > capacity(out(2))
                         fprintf(route, '\n %s \n', 'Il giro è formato dal nodo out(2):  ');                          
                         fprintf(route, '\n %s', 'Si può fare il giro con il truck di capacità:  '); 
                         fprintf(route, '%d \n', trucks(i)); 
                         ifatto=i;
                         break
                     end
                 end
              
                 for i=2:n_trucks 
                     if trucks(i) > capacity(out(1)) && ifatto~=i
                         fprintf(route, '\n %s \n', 'Il giro è formato dal nodo out(1)');                         
                         fprintf(route, '\n %s', 'Si può fare il giro con il truck di capacità:  '); 
                         fprintf(route, '%d \n', trucks(i)); 
                         i1fatto=i;
                         break
                     end
                 end
              
                    if ifatto==0
                      fprintf(route, '\n %s \n', 'Il nodo out(2) non può essere servito da nessun trucks!'); 
                    end
                    
                    if i1fatto==0
                     fprintf(route, '\n %s \n', 'Il nodo out(1) non può essere servito da nessun trucks!'); 
                    end
                 
                 else
                     
                  for i=2:n_trucks
                    if trucks(i) >= capacity(out(1))
                       fprintf(route, '\n %s \n', 'Il giro è formato dal nodo out(1)');                        
                       fprintf(route, '\n %s', 'Si può fare il giro con il truck di capacità:  '); 
                       fprintf(route, '%d \n', trucks(i)); 
                       ifatto=i;
                       break
                    end
                  end
                 
                  for i=2:n_trucks 
                       
                         t=i;
                         y=capacity(out(2));
                         z=trucks(i);
                         
                     if trucks(i) > capacity(out(2)) && ifatto~=i
                        fprintf(route, '\n %s \n', 'Il giro è formato dal nodo out(2)');            
                        fprintf(route, '\n %s', 'Si può fare il giro con il truck di capacità:  '); 
                        fprintf(route, '%d \n', trucks(i)); 
                        i1fatto=i;
                        break
                     end
                     
                  end
                     
                 end
                 
                 if ifatto==0
                      fprintf(route, '\n %s \n', 'Il nodo out(1) non può essere servito da nessun trucks!'); 
                 end
                 
                 if i1fatto==0
                     fprintf(route, '\n %s \n', 'Il nodo out(2) non può essere servito da nessun trucks!'); 
                 end
                     
                 end
   
                 
                 
                 
                 
     
                    
                     
                   
         fclose(route);
        end
        
         n_point=1 %se size(data,1)<1 &&  n_trucks=0 && n_point<1, significa che ho raggiunto tutti i nodi ed ho solo il DC
         
    end
    
end
end

function [tour, length,out] = clarke_wright(datain, capacity, truck,n, num)

out=[] %creo il vettore che conterrà i punti non raggiunti dal route che la function sta per creare

capacity_tot=truck   %considero la capacità del PRIMO mezzo

hubIdx = 1;     %il distribution center è il primo punto 

distances=datain; 

savings = zeros(n);  % creo una matrice dove scrivere la matrice di SAVINGS

for i=1:n  %calcolo la matrice di savings
    
    if i==hubIdx
       continue;
    end
    
    savings(i,(i+1):n)=distances(i,hubIdx)+distances(hubIdx,(i+1):n)-distances(i,(i+1):n);

end

sav=savings

minParent = 1:n;

[~,si] = sort(savings(:),'descend'); %ordino i savings ottenuti in maniera decrescente
si=si(1:ceil(end/2));                % posizione all'interno della matrice 

Vh = zeros(1,n);  %vettore che mi consente di vedere i punti che tolgo durante la costruzione del giro
Vh(hubIdx) = 1;
VhCount = n-1;   %serve per contare i punti ancora rimasti da inserire nel giro
degrees = zeros(1,n);  %vettore per i gradi di ogni nodo. il grado indica i collegamenti con altri nodi. ovviamente quando il grado è pari a 2, non posso più utilizzare quel nodo.

selectedIdx = 1;  % edge to try for insertion

tour = zeros(n,2);  %creo il tour indicando gli estremi dell'edge inserito
curEdgeCount = 1;
cap_disponibile=capacity_tot; %capacità iniziale (pari a quella del mezzo)

fg=0;   %numero nodi che rimangono fuori dal giro di consegna


while VhCount>2 %ho ancora nodi da inserire (VhCount = 2 signifca che ho il Dc e un solo nodo; il giro è quindi chiaro)
      
   %if si(selectedIdx) == 1   %significa che il saving è nullo e quindi dobbiamo uscire dal while!!!
      
   %individuo la colonna e la riga riferite al valore del saving
   %considerato. colonna e riga si riferiscono ai 2 nodi che permettono
   %tale saving.
   
    i = mod(si(selectedIdx)-1,n)+1       %colonna 
    j = floor((si(selectedIdx)-1)/n)+1   %riga

    % m=Vh  
  
    if Vh(i)==0 && Vh(j)==0 && (minParent(i)~=minParent(j)) && i~=j && i~=hubIdx && j~=hubIdx    % always all degrees are <= 2, so it is not required to check them
       %g=selectedIdx

       carico=0; %carico iniziale
       
       if degrees(i)==0 
          carico=capacity(i)  %aggiungo al carico, la capacità del primo punto
       end
            
       if degrees(j)==0
          carico= carico +capacity(j)  %aggiungo al carico, la capacità del secondo punto
       end
           
       %capdisp=cap_disponibile
         
       if cap_disponibile >= carico   %verifico che la quantità che carico non superi la capacità tot del mezzo
                            
          cap_disponibile= cap_disponibile - carico  %capacità residua 
         
          %incremento i gradi dei 2 punti che ho inserito
          
          degrees(i)=degrees(i)+1;
          degrees(j)=degrees(j)+1;
          
          %aggiungo i 2 nodi al tour
          
          tour(curEdgeCount,:) = [i,j];

          if minParent(i)<minParent(j)
             minParent(minParent==minParent(j))=minParent(i);
          else
             minParent(minParent==minParent(i))=minParent(j);            
          end

 
          curEdgeCount = curEdgeCount + 1; 
            
            
          if degrees(i)==2   %il nodo non può più essere inserito
             Vh(i) = 1; 
             VhCount = VhCount - 1;  %tolgo il nodo dai nodi ancora da inserire
          end

          if degrees(j)==2   %il nodo non può più essere inserito
             Vh(j) = 1;
             VhCount = VhCount - 1;  %tolgo il nodo dai nodi ancora da inserire
          end
              
              
              
     else if degrees(i) == 0  %se il carico dei nodi i + j è maggiore della capacità disponibile, 
                              %verifico quale dei due nodi posso inserire.
                              %se j è già stato inserito in precedenza
                              %(degrees(j)=1),vado a togliere il nodo i che
                              %non è mai stato inserito e lo scrivo nel
                              %vettore [out].
                  
          if degrees(j) == 1
             Vh(i) = 1;  %tolgo i
             VhCount = VhCount - 1;
       
             fg=fg+1;  %incremento i nodi fuori dal giro
             %bugi = input('punto i  ')
             %bugi=i
               
             out= [out ;i];  %inserisco nel vettore [out] il punto i che rimane fuori
       
         else if cap_disponibile<capacity(i)  % tolgo il nodo i poichè la sua capacità è > di quella disponibile
            
                 Vh(i) = 1;   %tolgo i
                 VhCount = VhCount - 1;
       
                 fg=fg+1;  %incremento i nodi fuori dal giro
                 % bugi = input('punto i  ');
                 %bugi=i;
               
                 out= [out ;i];  %inserisco nel vettore [out] il nodo i che rimane fuori
             
             else               %la capacità di i < capacità disponibile. posso inserire il nodo i e devo togliere j.
                 Vh(j) = 1;     %tolgo j
                 VhCount = VhCount - 1;
                 fg=fg+1; %incremento i punti fuori dal giro
                 %bugj = input('punto j  ');
                 %bugj=j;
                 out= [out ;j];   %inserisco nel vettore [out] il nodo j che rimane fuori
            
             end
          end
          
         else            % il nodo i ha degrees = 1. quindi vado ad inserire tale nodo i e tolgo il nodo j.
            Vh(j) = 1;     %tolgo j
            VhCount = VhCount - 1;
            fg=fg+1; %incremento i punti fuori dal giro
            %bugj = input('punto j  ');
            %bugj=j;
            out= [out ;j];   %inserisco nel vettore [out] il nodo j che rimane fuori
            
        end

                   
     end
             
  end
    
selectedIdx = selectedIdx + 1
   %else
   % vh = 2   verificare la capacità. altrimenti Vh == 0 e fa fuori giro. 
end
      
%inserisco i punti estremali del giro: il primo punto dopo la partenza dal
%DC e l'ultimo punto prima di ritornare al DC

remain = find(Vh==0);  %individuo i nodi non inseriti
n1=remain(1)
n2=remain(2)

  if fg>0
     tour(n-fg+1:n,:)= []  %il tour comprende i nodi totali meno quelli fuori giro + il DC.
  end
  
  v=degrees(:,:);
  
if isequal (degrees, zeros (1,n) )  %se c'è uguaglianza: non ha inserito nemmeno un punto! devo verificare che remain 1 e 2 non superino la capacità disponibile del mezzo.
 
    if capacity(n1) + capacity (n2) <= cap_disponibile    %posso inserire questo arco
       giro= zeros (1,3); 
       giro (1) = 1;  %parto dal DC
    if n1 > n2   
       giro (2) = n2; 
       giro (3) = n1;
    else
       giro (2) = n1;
       giro (3) = n2;
    end
    
    tour = giro;  %il tour è formato dal DC e i 2 nodi rimasti
    length = distances(n1,n2) + distances(1,n2) + distances(1,n1); %lunghezza del tour
   
    else if capacity (n2) < capacity (n1) %inserisco nel giro solo un nodo, dando precedenza a quello con capacità maggiore.
            tour = [1 n1];
            out =  [out ;n2];
            length= 2 * distances(1,n1);
        else 
            tour = [1 n2];
            out =  [out ;n1];
            length= 2 * distances(1,n2);
        end
    end
else  
    
tour(curEdgeCount,:) = [hubIdx n1];  %primo punto 
curEdgeCount = curEdgeCount + 1;
tour(curEdgeCount,:) = [hubIdx n2]; %ultimo punto

tour = stitchTour(tour)  %uso la function per ordinare i punti costituenti il giro
tour=tour(:,1)';
length=distances(tour(end),tour(1));  %distanza nodo finale - DC

  for i=1:n-(1+fg)  
      length=length+distances(tour(i),tour(i+1));  %somma delle distanze tra i nodi
  end

end
  

%scrivo su file i giri e i punti che rimangono out e i successivi giri che
%si creano con gli out dei giro precedente

route=fopen('route10.txt', 'at');
fprintf(route, '\n %s', 'I nodi del route sono raggiungi dal mezzo di capacità:  ');
fprintf(route, '%d \n', truck);
fprintf(route, '\n %s \n', 'Il route creato è formato dai nodi. NB:dopo l''ultimo nodo si ritorna al nodo 1. [route]');
fprintf(route, '%d \n', tour);
fprintf(route, '\n %s \n', 'Corrispondenti ai nodi iniziali di giro');
fprintf(route, '%d \n', num(tour));

fprintf(route, '%s \n', 'Nodi non inseriti [out]');
fprintf(route, '%d \n', out); 
fprintf(route, '\n %s \n', 'Corrispondenti ai nodi iniziali di giro');
fprintf(route, '%d \n', num(out));



fclose(route);



end



function tour = stitchTour(t) % uniforma il tour [a b; b c; c d; d e;.... ]

n=size(t,1);

[~,nIdx] = sort(t(:,1));
t=t(nIdx,:);

tour(1,:) = t(1,:);
t(1,:) = -t(1,:);
lastNodeIdx = tour(1,2);

for i=2:n
    nextEdgeIdx = find(t(:,1)==lastNodeIdx,1);
    if ~isempty(nextEdgeIdx)
        tour(i,:) = t(nextEdgeIdx,:);
        t(nextEdgeIdx,:)=-t(nextEdgeIdx,:);
    else
        nextEdgeIdx = find(t(:,2)==lastNodeIdx,1);
        tour(i,:) = t(nextEdgeIdx,[2 1]);
        t(nextEdgeIdx,:)=-t(nextEdgeIdx,:);
    end
    lastNodeIdx = tour(i,2);
end


end