N=51;                                                                       %Arithmos pelaton
Qmax=160;                                                                   %Xoritikotita oximatos
Tmax=200;                                                                   %Megisto mikos diadromis
XE=10;                                                                      %Xronos eksipiretisis
Nxy=[30 40;37 52;49 49;52 64;20 26;40 30;21 47;17 63;31 62;52 33;51 21;42 41;31 32;5 25;12 42
    36 16;52 41;27 23;17 33;13 13;57 58;62 42;42 57 ; 16 57 ; 8 52 ; 7 38 ; 27 68;30 48;43 67
    58 48;58 27;37 69;38 46;46 10;61 33;62 63;63 69;32 22;45 35;59 15;5 6;10 17;21 10;5 64
    30 15;39 10;32 39;25 32;25 55;48 28;56 37];                             %Syntetagmenes/Pinakas kostoys euklidies apostaseis

di=[0;7;30;16;9;21;15;19;23;11;5;19;29;23;21;10;15;3;41;9;28;8;8;16;10;28;7;15;14;6;19;11;12;23
    26;17;6;9;15;14;7;27;13;11;16;10;5;25;17;18;10];                        %zitisi

for i=1:N                                                                   %Dimioyrgia apostasewn
    for j=1:N
        if(i==j)
            C(i,j)=0;                                                       %H apostasi toy kathe simeioy apo ton eayto toy einai mhden
        else
            C(i,j)=((Nxy(i,1)-Nxy(j,1))^2+(Nxy(i,2)-Nxy(j,2))^2)^(1/2);     %Pinakas apostaseon
        end
    end
end
standarC=zeros(N,N);                                                        %pinakes bohthitikoi pou kratane ton arxiko pinaka apostasewn
extrastandarC=zeros(N,N);
for i=1:N
    for j=1:N
        standarC(i,j)=C(i,j);
        extrastandarC=C(i,j);
    end
end

MAXDSTNC=-1;                                                                %boithitiki metabliti
for i=1:N
    for j=1:N
        if (standarC(i,j)>MAXDSTNC)
            MAXDSTNC=standarC(i,j);                                         %boithitiki metabliti (h megalyterh apostasi ston pinaka)
        end
    end
end

X=zeros(N,N);                                                               %Mhdenismos pinaka poy mas deixnei poioys eksyphrethsa
Qol=0;                                                                      %Synoliko Qol gia kathe diadromi (xwthtikothta oximatos)
xrol=0;                                                                     %Synoliko xrol gia kathe diadromi (xronos paramonis oximatos sthn diadromi)
Q=0;                                                                        %Ypologismos prosorinoy Q gia kathe mia diadromi
xr=0;                                                                       %Ypologismos prosorinoy xr gia kathe diadromi
count=1;                                                                    %arithmos diadromwn (monopatiwn)
Xdiadromes=zeros(N,N);                                                      %Mhdenismos pinaka poy mas deixnei poioys eksyphrethsa ana diadromi (monopatia)
cost=0;                                                                     %beltisto kostos/trexon kostos
sum=0;
position=1;                                                                 %boithikikh metavliti pou mou deixnei se poion pelati vriskomai (1=apothiki)
clientdistance=0;
client=0;
for n=1:N
    for m=1:N
        if(C(n,m)>clientdistance)
            clientdistance=C(n,m);
        end
    end
end
maxclientdistance=clientdistance+10;
maxdstnc=maxclientdistance;
while (sum<N-1)                                                             %VRP (Algorithmos plhsiesterou geitona me periorismous xwrhtikothtas kai paramonis oximatos sthn diadromi(62-107))
    sum=0;                                                                  %metavliti pou deixnei kathe stigmi posoi pelates exoun eksipiretithei (termatiki sinthiki loupas)
    maxclientdistance=maxdstnc; 
    for j=2:N                                                               %bronxos poy entopizetai o pio kodinos pelaths ap to position (65-71)
        if((j~=position)&&(C(position,j)~=-1)&&(C(position,j)<maxclientdistance))               
            maxclientdistance=C(position,j);
            clientdistance=C(position,j);
            client=j;
        end
    end
    
    if ((Q+di(client)<=Qmax) && (xr+XE+C(position,client)+C(client,1)<=Tmax))   %elenxw periorismous wste na dw an mporw na kanw thn kinhsh pros ton client pou entopisa prohgoumenws (sto xr perilamvanetai kai o xronos epistrofis stin apothiki apo to shmeio pou vriskomai) (73-92)
        xr=xr+XE+C(position,client);                                        %an mporw tote ypologizw to neo xr kai to neo Q  
        Q=Q+di(client);
        X(position,client)=1;                                               %shmeiwnw me asso son arxika zeroPinaka to shmeio poy ypodeilwnei thn kinhsh ap to position ston client
        Xdiadromes(position,client)=count;                                  %antistoixa shmeiwnw me ton arithmo ths yparxousas diadromis ton antistoixo pinaka twn diadromwn
        for i=1:N                                                           %tis metakinhseis pleon pros kai apo ton pelati pou molis eksyphretithike tis diagrafw
            C(i,client)=-1;
        end
        position=client;
    else                                                                    %an den borw na kanw tin kinisi 
        xr=xr+C(position,1);                                                %enhmerwnw to xr mou oti gyrnaei stin apothiki
        Qol=[Qol;Q];                                                        %ekxwrw to q ths oloklirwmenis diadromis
        xrol=[xrol;xr];                                                     %ekxwrw to xr ths oloklirwmenis diadromis
        X(position,1)=1;
        Xdiadromes(position,1)=count;                                       %enhmwrwnw tis kinhseis stous pinakes X ap to client sthn apothiki
        xr=0;                                                               %mhdenizw xwritikotites kai xronous diadromis gia epomenh xrhsh
        Q=0;
        count=count+1;                                                      %enhmerwnw gia thn dimiourgia tis epomenis diadromis
        position=1;                                                         %ksana orizw ws shmeio ekkinisis thn spothiki
    end
    
    for n=1:N                                                               %ypologizw se kathe epanalipsi to sum wste an exoun eksypiretithei 50 pelates na oloklirwthei to provlima
        for m=2:N
            sum=sum+X(n,m);
        end
    end
    if (sum==N-1)                                                           %termatikos elegxos oloklirwsis tis diadikasias kai ekxwrhshw twn telautaiwn q,xr
         xr=xr+C(position,1);
         Qol=[Qol;Q];                               %Apothikeyoyme ta temaxia poy edose to fortigo
         xrol=[xrol;xr];                            %Apothikeyoyme to mhkos diadromhs
         X(position,1)=1;                                  %Gyrname sthn apothiki
         Xdiadromes(position,1)=count;
    end
end

for i=1:N                                                                   %YPOLOGISMOS KOSTOUS VRP
    for j=1:N
        cost=cost+standarC(i,j)*X(i,j);                                     
    end
end                                                                         

cost

Xcordinates= Nxy(1,1);
Ycordinates= Nxy(1,2);
max=count;
startpoint=-1;
for i=1:count
    startpoint=1;
    acefound=1;
    while (acefound~=2)
        for j=1:N
            if (Xdiadromes(startpoint,j)==i)
          
                Xcordinates= [Xcordinates;Nxy(j,1)];
                Ycordinates= [Ycordinates;Nxy(j,2)];
                startpoint=j;
                j=0;
                if (startpoint==1)
                    if(i==1) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-r')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==2) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-g')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==3) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-b')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==4) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-y')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==5) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-m')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==6) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-c')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==7) && (i<=max)
                        plot(Xcordinates,Ycordinates,'o-k')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i>=8)
                        plot(Xcordinates,Ycordinates,'o-r')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    
                end
            end
        end
                        if (i==max)
                            hold off
                            acefound=2;
                            break
                        end
    end
end

X2diadromes=Xdiadromes;                                                     %Arxikopoihsh boithitikwn metavlitwn gia ta relocate,exchange
X2=X;
repetitions=0;
a=0;
aproig=0;
aepom=0;
b=0;
bepom=0;
firstroute=0;
secondroute=0;
Q1=0;
Q2=0;
xr1=0;
xr2=0;
cost2=0;

tabuvar1=0;
bttrcstfnd=0;                                                               % h metabliti bttrcstfnd (better cost found) metraei poses fores exei beltiwthei to kostos ap to arxiko

while(bttrcstfnd~=15)                                                       %RELOCATE 1-0 (135-217) algorithmos topikis anazitisis relocate 1-0 
    X2diadromes=Xdiadromes;                                                 %X2 dokimastikoi pinakes diadromwn pou ypokeintai se allagi logw tis topikis anazitisis
    X2=X;
    
    a=randi([2,N]);                                                        %dinoume tyxaia timi gia enan pelati gia to relocate 1-0
    for j=1:N                                                              %vriskoume ton proigoumeno tou                             
        if (X2(j,a)==1)            
            aproig=j;             
        end
    end
            
    for j=1:N                                                              %vriskoume ton epomeno tou
        if (X2(a,j)==1)             
            aepom=j;                
        end
    end
    firstroute=X2diadromes(aproig,a);                                       %vriskoume ton arithmo tou monopatiou sto opoio anhkei
    
    b=randi([2,N]);                                                        %dinoume tyxaia timi se enan allon pelati gia to relocate
    while(b==a || b==aproig)                                                %den theloume o pelaths b na einai o idios me ton a
        b=randi([2,N]);
    end
    for j=1:N
        if (X2(b,j)==1)                                                     %vriskoume ton epomeno tou b
            bepom=j;           
        end
    end
    secondroute=X2diadromes(b,bepom);                                       %vriskoume ton arithmo tou monopatiou tou pelati b
    
    Q1=0;                                                                   
    Q2=0;
    xr1=0;
    xr2=0;
    
    X3diadromes=X2diadromes;                                                %allazoume stous pinakes tis sindeseis metaksi pelatwn a kai b
    X2diadromes(aproig,a)=0;              
    X2diadromes(a,aepom)=0;
    X2diadromes(aproig,aepom)=firstroute;
    X2diadromes(b,bepom)=0;
    X2diadromes(b,a)=secondroute;
    X2diadromes(a,bepom)=secondroute;
    
    for i=1:N                                                               %vriskoume ta nea q kai xr gia tia nea monopatia pou dimiourgithikan (177-188)
        for j=1:N
            if (X2diadromes(i,j)== firstroute)   
                Q1=Q1+di(j);                   
                xr1=xr1+XE+standarC(i,j);            
            end
            if (X2diadromes(i,j)== secondroute)  
                Q2=Q2+di(j);                    
                xr2=xr2+XE+standarC(i,j);            
            end
        end
    end
    
    if (Q1<=Qmax && xr2<=Tmax && Q2<=Qmax && xr2<=Tmax)                     %elegxw an ta nea q kai xr ikanopoioun tous periorismous
        X3=X2;                                                              %monimopoiw tin allagi an auth exei ginei apodekti
        X2(aproig,a)=0;              
        X2(a,aepom)=0;
        X2(aproig,aepom)=1;
        X2(b,bepom)=0;
        X2(b,a)=1;
        X2(a,bepom)=1;
        
        for i=1:N                                                          %ypologizw to kostos tis neas diadromis cost2
             for j=1:N
                cost2=cost2+standarC(i,j)*X2(i,j);   
             end
        end
        if (cost2<cost)                                                    %elegxw an to neo kostos einai kalytero toy prohgoumenou
            X=X2;                                
            Xdiadromes=X2diadromes;
            cost=cost2;
            bttrcstfnd=bttrcstfnd+1;
        end
        cost2=0;
    else
        X2diadromes=X3diadromes;                                            %
    end
    
    
    
end

Xcordinates= Nxy(1,1);
Ycordinates= Nxy(1,2);
max=count;
startpoint=-1;
for i=1:count
    startpoint=1;
    acefound=1;
    while (acefound~=2)
        for j=1:N
            if (Xdiadromes(startpoint,j)==i)
          
                Xcordinates= [Xcordinates;Nxy(j,1)];
                Ycordinates= [Ycordinates;Nxy(j,2)];
                startpoint=j;
            
                if (startpoint==1)
                    if(i==1) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-r')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==2) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-g')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==3) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-b')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==4) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-y')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==5) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-m')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==6) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-c')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==7) && (i<=max)
                        plot(Xcordinates,Ycordinates,'o-k')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i>=8)
                        plot(Xcordinates,Ycordinates,'o-r')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    
                end
            end
        end
                        if (i==max)
                            hold off
                            acefound=2;
                            break
                        end
    end
end

cost

bttrcstfnd1=0;                                                              %arxikopoihseis metavlitwn gia to exchange paromoioi me autous tou relocate1-0 (221-231)
extndcost=0;
w=0;                                                                        
wproig=0;                                                                   
wepom=0;                                                                    
l=0;                                                                        
lproig=0;                                                                   
lepom=0;                                                                    
routeone=0;
routetwo=0;
cost3=0;

LIST=zeros(11,2);                                                           %arxikopoihsh metavlitwn gia to tabuList (233-241)
LSTSIZE=10;                                                                 %megethos tabuList
currentsize=0;
tabuvar1=0;                                                                 %tabuvar1 deixnei an h lysh pou eksetazetai einai kalyterh
tabuvar2=0;                                                                 %tabuvar2 deixnei an h eksetazomenh kinhsh yparxei h oxi sthn lista
tabuvar3=0;                                                                 %tabuvar3 deixnei an h eksetazomenh lyei einai kalyterh ths twrinhs dexomenoi ena 10 tois ekato apoklisi
tabuvar2temp=0;
tabuvar3rptfail=0;                                                          %deixnei poses epanalipseis ginane pou dwsane mi epithimito apotelesma
tabuvar4rptpos=0;                                                           %deixnei poses epanalipseis ginane apodektes

Xdiag=X;                                                                    %pinakes Xdiag aforoun tis diadikasies ths diaforopoihshs ka entatikopoihshs kai einai oi X pinakes mou se anw trigwnikous logw symmetrikotitas
Xdiaghlp=zeros(N,N);
total=0;
for i=1:N
    j=1;
    while(i<j && j<=N)
        Xdiag(i,j)=Xdiag(i,j) + Xdiag(j,i);
        j=j+1;
    end
end
for i=1:N
    j=1;
    while(i>j)
        Xdiag(i,j)=0;
        j=j+1;
    end
end

while(tabuvar4rptpos~=20000)                                                %o brogxos autos emperiexei thn topiki anaziti exchange ka8ws kai ton TabuSearch (261-728)
    X2diadromes=Xdiadromes;                                                 %o Exchange 1-1 algorithmos ginetai stis 262-373
    X2=X;
    
    w=randi([2,N]);                                                         %paromoia diadikasia me auth tou Relocate 1-0
    for j=1:N
        if (X2(j,w)==1)                                                    
            wproig=j;                                                     
        end
    end
    for j=1:N
        if (X2(w,j)==1)                                                    
            wepom=j;                                                        
        end
    end
    routeone=X2diadromes(wproig,w);                                         
    
    l=randi([2,N]);                                                       
    while(l==w)                                                             
        l=randi([2,N]);
    end
    for j=1:N
        if (X2(j,l)==1)                                                     
            lproig=j;                                                    
        end
    end
    for j=1:N
        if (X2(l,j)==1)                                                     
            lepom=j;                                                       
        end
    end
    routetwo=X2diadromes(lproig,l);
    
    Q1=0;
    Q2=0;
    xr1=0;
    xr2=0;
    
    if (l==wproig)                                                          %paromoia diadikasia me ton Relocate 1-0 me thn methodologia toy Exchange 1-1
        X3diadromes=X2diadromes;
        X2diadromes(wproig,w)=0;                                           
        X2diadromes(w,wepom)=0;
        X2diadromes(lproig,l)=0;
        X2diadromes(lproig,w)=routeone;
        X2diadromes(w,l)=routeone;
        X2diadromes(l,wepom)=routeone;          
    elseif(w==lproig)
        X3diadromes=X2diadromes;     
        X2diadromes(wproig,w)=0;                                           
        X2diadromes(w,l)=0;
        X2diadromes(l,lepom)=0;
        X2diadromes(wproig,l)=routeone;
        X2diadromes(l,w)=routeone;
        X2diadromes(w,lepom)=routeone;
    else 
        X3diadromes=X2diadromes;          
        X2diadromes(wproig,w)=0;                                           
        X2diadromes(w,wepom)=0;
        X2diadromes(lproig,l)=0;
        X2diadromes(l,lepom)=0;
        X2diadromes(wproig,l)=routeone;
        X2diadromes(l,wepom)=routeone;
        X2diadromes(lproig,w)=routetwo;
        X2diadromes(w,lepom)=routetwo;
    end
    
    for i=1:N                                                               %ypologismos q kai xr gia tis nees diadromes
        for j=1:N
            if (X2diadromes(i,j)== routeone)   
                Q1=Q1+di(j);                   
                xr1=xr1+XE+standarC(i,j);            
            end
            if (X2diadromes(i,j)== routetwo)  
                Q2=Q2+di(j);                    
                xr2=xr2+XE+standarC(i,j);            
            end
        end
    end
    
    if (Q1<=Qmax && xr1<=Tmax && Q2<=Qmax && xr2<=Tmax)                     %an oi periorismoi ikanopoiountai 
        if(l==wproig)
            X3=X2;
            X2(wproig,w)=0;                                                
            X2(w,wepom)=0;
            X2(lproig,l)=0;
            X2(lproig,w)=1;
            X2(w,l)=1;
            X2(l,wepom)=1;          
        elseif(w==lproig)
            X3=X2;     
            X2(wproig,w)=0;                                                 
            X2(w,l)=0;
            X2(l,lepom)=0;
            X2(wproig,l)=1;
            X2(l,w)=1;
            X2(w,lepom)=1;
        else 
            X3=X2;          
            X2(wproig,w)=0;                                                 
            X2(w,wepom)=0;
            X2(lproig,l)=0;
            X2(l,lepom)=0;
            X2(wproig,l)=1;
            X2(l,wepom)=1;
            X2(lproig,w)=1;
            X2(w,lepom)=1;
        end
         
         for i=1:N                                                          %ypologismos kostous lyshs tou exchange 1-1                                     
             for j=1:N
                cost3=cost3+standarC(i,j)*X2(i,j);   
             end
         end
         
         if (cost3<cost)                                                    %an to kostos einai kalytero tou twrinou to ekxwrw ws to neo beltisto
            X=X2;                                                           
            Xdiadromes=X2diadromes;
            cost=cost3;
            cost
            tabuvar1=1;                                                     %ekxwrw ston metriti oti vrethike kostos kalytero                                                    
            bttrcstfnd1=bttrcstfnd1+1;
         end 
         
         if(cost3<(cost+(0.1*cost)))                                        %elegxw an to kostos ths lyshs moy einai kalytero apo thn mexri twra veltisti exwdas prosuesei se ayth ena epipleon deka tois ekto tis timis tis
             tabuvar3=1;
             extndcost=extndcost+1;
         end
         
        cost3=0;
                                                                            %ekkinhsh diadikasias TABUSEARCH 391-728
        for i=1:10                                                          %TABU prosedure starts,elegxw an h kinhsh moy yparxei sthn lista
            if((LIST(i,1)==w && LIST(i,2)==l) || (LIST(i,1)==l && LIST(i,2)==w))     %energopoiw thn metavliti tabuvar2 an h kinhsh yparxei
                tabuvar2temp=i;
                tabuvar2=1;                                                
                break
            end
        end
        
        if(tabuvar2==0)                                                     %an h kinhsh mou den yparxei thn ekxwrw                                   
            for i=10:-1:1 
                if (currentsize==0)
                    currentsize=currentsize+1;
                    break
                end
                currentsize=currentsize+1;
                for j=1:2
                    LIST(i+1,j)=LIST(i,j);
                end
            end
            LIST(11,1)=0;
            LIST(11,2)=0;
            LIST(1,1)=w;
            LIST(1,2)=l;
        else                                                                %alliws ean yparxei elegxw an einai kalyterh h oxi ths prohgoumenhs
            if(tabuvar1==1)                                                 %an h lysh einai kalyterh ths proigoumenis
                push=push+1;
                if(tabuvar2temp<3)
                    if(tabuvar2temp==2)                                     %thn krataw kai thn ekxwrw kanonika
                        LIST(2,1)=LIST(1,1);
                        LIST(2,2)=LIST(1,2);
                        LIST(1,1)=w;
                        LIST(1,2)=l;
                    end
                else
                    for i=(tabuvar2temp-1):-1:1
                        for j=1:2
                            LIST(i+1,j)=LIST(i,j); 
                        end
                    end
                    LIST(1,1)=w;
                    LIST(1,2)=l;
                end
            else                                                            %an yparxei kai ean h lysh den einai kalyterh ths proigoymenis den ephreazetai to programma kai apla aporriptetai
                    
                 
            end
        end
        
        if(tabuvar3==1)                                                     %(!!!!)an h lysh mou einai kalyterh apo thn beltisti estw kai logw ths ayjhshs ths kata deka tois ekato ekxwreite h kinhsh ston trigvniko pinaka diadromwn poy ua xrhsimopoijuei gia tin diaforopoihsh kai thn entatikopoihsh
            for i=1:N
                for j=1:N
                    Xdiag(i,j)=Xdiag(i,j)+X2(j,i)+X2(i,j);
                    if (i>=j)
                        Xdiag(i,j)=0;
                    end
                end
            end
        end
        
        tabuvar2temp=0;
        tabuvar4rptpos=tabuvar4rptpos+1;                                    %auksanw thn metavliti tabuvar4rptpos se kathe lysi poy einai mesa sto orio kai katw tou 10% ths veltistis
    else
        tabuvar3rptfail=tabuvar3rptfail+1;
        X2diadromes=X3diadromes;
    end
    tabuar3=0;                                                              %mhdenismos twn arxikwn metavlitwn tou tabusearch
    tabuvar2=0;
    tabuvar1=0;
    
    if((tabuvar4rptpos==150) || (tabuvar4rptpos==300) || (tabuvar4rptpos==450) ||(tabuvar4rptpos==600)||(tabuvar4rptpos==750)||(tabuvar4rptpos==1000))      % (!!!)Brogxos ths diadikasias edadikopoihshs 460-574
        if((tabuvar4rptpos==150) || (tabuvar4rptpos==300) || (tabuvar4rptpos==450) ||(tabuvar4rptpos==600)||(tabuvar4rptpos==750)||(tabuvar4rptpos==1000))
            tabuvar4rptpos=tabuvar4rptpos+1;
        end                                                                 %h if ayth epitrepei thn diadikasia tis edadikopoihshs na symbainei mia fora ana 150 apodektes lyseis tou provlimatos
        entatic=Xdiag;                                                      %arxikopoihsh metavlitwn pinakwn X gia thn edatikopoihsh
        Xentatic=zeros(N,N);
        Xentaticdiadromes=zeros(N,N);
        cost4=0;                                                            %kostos pou borei na prokypsei apo edatikopoihsh
        
        grammi=1;
        maxclient=1;
        maxroutes=0;
        Qa=0;
        xrA=0;
        xrB=0;
        routenumber=1;
        previousmaxclient=0;
        total=0;
        while(total~=N-1)                                                   %mexri na syblirwthei to synolo twn 50pelatwn
            maxroutes=0;
            if (grammi~=N)
                for j=grammi:N                                              %vriskw ton pio dimofili pelath apo position 'grammi' (480-499)
                    if (j==grammi)
                        maxclient=j;
                        maxroutes=entatic(grammi,j);
                    end
                    if(entatic(grammi,j)>maxroutes)
                        maxclient=j;
                        maxroutes=entatic(grammi,j);
                    end
                end
            end
            if(grammi~=1)
                for i=1:grammi
                    if(entatic(i,grammi)>maxroutes)
                        maxclient=i;
                        maxroutes=entatic(i,grammi);
                    end
                end
            end
            
            if(maxclient==1)                                                %se periptwsi pou o dimofiloteros pelaths einai h apothiki toy systhmatos mou
                Xentaticdiadromes(grammi,maxclient)=routenumber;            %katagrafetai h kinhsh stous pinakes X
                Xentatic(grammi,maxclient)=1;
                if(grammi~=1)                                               %diagrafodai oi kinhseis toy  prwhn dimofilesterou pelath me kathe allon
                    for i=grammi:N
                        entatic(grammi,i)=0;
                    end
                    for i=1:grammi
                        entatic(i,grammi)=0;
                    end
                end
                routenumber=routenumber+1;                                  %kai enhmerwnondai oi metavlites q,xr kai arithmou monopatiou routenoumber
                grammi=1;
                Qa=0;
                xrA=0;
                xrB=0;
            else                                                            %an o dimofilesteros den einai h apothiki
                Qa=Qa+di(maxclient);                                        %enhmerwnw tis metavlites q kai xr twn diadromwn pou ephreazodai
                xrA=xrA+XE+standarC(grammi,maxclient);
                xrB=xrA+standarC(maxclient,1);
                if((Qa<=Qmax) && (xrB<=Tmax))                               %kai an oi metavlites ikanopoioun tous periorismous katagrafw thn kinhsh
                    Xentaticdiadromes(grammi,maxclient)=routenumber;
                    Xentatic(grammi,maxclient)=1;
                    total=total+1;
                    if(total==(N-1))
                        Xentaticdiadromes(maxclient,1)=routenumber;
                        Xentatic(maxclient,1)=1;
                    end
                    if(grammi~=1)                                           %diagrafodai oi syndeseis toy pelath poy eksyphreteitai me opoiondipote allo pelath-komvo
                        for i=grammi:N
                            entatic(grammi,i)=0;
                        end
                        for i=1:grammi
                            entatic(i,grammi)=0;
                        end
                    end
                    entatic(grammi,maxclient)=0;                            %kai pleon o neos pelaths grammi ginetai aytos poy molis ekgyphreththhke
                    grammi=maxclient;

                else                                                        %an oi periorismoi den ikanopoioudai tote to oxhma mou epistrefei sthn apothiki
                    Xentaticdiadromes(grammi,1)=routenumber;
                    Xentatic(grammi,1)=1;
                    routenumber=routenumber+1;
                    if(grammi~=1)
                        for i=grammi:N
                            entatic(grammi,i)=0;
                        end
                        for i=1:grammi
                            entatic(i,grammi)=0;
                        end
                    end
                    grammi=1;                                               %kai arxikopoiw ksana tis metavlites ths diadromhs mou
                    Qa=0;
                    xrA=0;
                    xrB=0;
                end
            end

        end
         for i=1:N                                                          %molis eksyphreththei olo to synolo twn pelatvn mou ypologizw to cost4 poy exei prokipsei apo thn edatikopoihsh                                                        
             for j=1:N
                cost4=cost4+standarC(i,j)*Xentatic(i,j);   
             end
         end
         cost4
         
         if (cost4<cost)                                                    %kai an to kostos einai kalytero tou veltistou to krataw ws neo beltisto
            X=Xentatic;                                                           
            Xdiadromes=Xentaticdiadromes;
            cost=cost4;                                                     
         end
         
         
    end
    
    if((tabuvar4rptpos==200) || (tabuvar4rptpos==400) || (tabuvar4rptpos==600) ||(tabuvar4rptpos==800)||(tabuvar4rptpos==1000)||(tabuvar4rptpos==1200))      % me thn if ayth ksekinaei h diadikasia ths diaforopoihshs
        if(tabuvar4rptpos==200) || (tabuvar4rptpos==400) || (tabuvar4rptpos==600) ||(tabuvar4rptpos==800)||(tabuvar4rptpos==1000)||(tabuvar4rptpos==1200)      
            tabuvar4rptpos=tabuvar4rptpos+1;
        end                                                                 %h diadikasia ths diaforopoihshs lemvanei xwra kath 200 dektes lyseis
        extrastandarC=zeros(N,N);
        for i=1:N
            for j=1:N
                extrastandarC=standarC(i,j);
            end
        end
        diafor=Xdiag;                                                       %arxikopoihsh metavlitwn gia thn diaforopoihsh
        Xdiafor=zeros(N,N);
        Xdiafordiadromes=zeros(N,N);
        cost5=0;
        seira=1;
        minclient=1;
        minroutes=1550;
        Qa=0;
        xrA=0;
        xrB=0;
        routenumber=1;
        total=0;
        for i=1:N                                                           %vgazoume ap tous pinakes tous mhdenikous arithmous epeidh ua anazitisoyme ton ligotero dimofili
            for j=1:N
                if (i==j)
                    diafor(i,j)=minroutes;
                end
            end
        end
        maxcldstnc=MAXDSTNC+1;
        maxst=MAXDSTNC+2;
        for i=1:N
            for j=1:N
                if (i==j)
                    extrastandarC(i,j)=maxst;
                end
            end
        end
        while (total~=(N-1))                                                %mexri na syblhrwthei to synolo twn pelatwn
            minroutes=1550-1;
            maxcldstnc=MAXDSTNC+1;
            if (seira~=N)                                                   %vriskoume ton ligotero dimofili (617-644)
                for j=seira:N
                    if (j==seira)
                        minclient=j;
                        minroutes=1550-1;
                        maxcldstnc=MAXDSTNC+1;
                    end
                    if((diafor(seira,j)<=minroutes) && (extrastandarC(seira,j)<=maxcldstnc))
                        minclient=j;
                        minroutes=diafor(seira,j);
                        maxcldstnc=extrastandarC(seira,j);
                    end
                end
            end
            if(seira~=1)
                for i=1:seira
                    if((seira==N) && (i==1))
                        minclient=i;
                        minroutes=1550-1;
                        maxcldstnc=MAXDSTNC+1;
                    end
                    if(diafor(i,seira)<=minroutes && (extrastandarC(i,seira)<=maxcldstnc))
                        minclient=i;
                        minroutes=diafor(i,seira);
                        maxcldstnc=extrastandarC(i,seira);
                    end
                end
            end
            
            if(minclient==1)                                                %an o ligotero dimofilis einai h apothiki
                Xdiafordiadromes(seira,minclient)=routenumber;
                Xdiafor(seira,minclient)=1;
                if(seira~=1)
                    for i=seira:N
                        diafor(seira,i)=0;
                    end
                    for i=1:seira
                        diafor(i,seira)=0;
                    end
                end
                routenumber=routenumber+1;
                seira=1;
                Qa=0;
                xrA=0;
                xrB=0;
            else                                                            %alliws an einai kanonikos pelaths
                Qa=Qa+di(minclient);                                        %enhmerwnoume to systhma me ta xarakthristika tou
                xrA=xrA+XE+standarC(seira,minclient);
                xrB=xrA+standarC(minclient,1);
                if((Qa<=Qmax) && (xrB<=Tmax))                               %an ikanopoiei tous periorismous katagrafetai sthn diadromi
                    Xdiafordiadromes(seira,minclient)=routenumber;
                    Xdiafor(seira,minclient)=1;
                    total=total+1;
                    if(total==(N-1))
                        Xdiafordiadromes(minclient,1)=routenumber;
                        Xdiafor(minclient,1)=1;
                    end
                    if(seira~=1)
                        for i=seira:N
                            diafor(seira,i)=1550;
                            extrastandarC(seira,i)=maxst;
                        end
                        for i=1:seira
                            diafor(i,seira)=1550;
                            extrastandarC(i,seira)=maxst;
                        end
                    end
                    diafor(seira,minclient)=1550;
                    extrastandarC(seira,minclient)=maxst;
                    seira=minclient;

                else                                                        %alliws to oxhma epistrefei sthn apothiki kai enhmerwnw pali tis adistoixes metavlites
                    Xdiafordiadromes(seira,1)=routenumber;
                    Xdiafor(seira,1)=1;
                    routenumber=routenumber+1;
                    if(seira~=1)
                        for i=seira:N
                            diafor(seira,i)=0;
                            extrastandarC(seira,i)=maxst;
                        end
                        for i=1:seira
                            diafor(i,seira)=0;
                            extrastandarC(i,seira)=maxst;
                        end
                    end
                    seira=1;
                    Qa=0;
                    xrA=0;
                    xrB=0;
                end
            end
            
        end
        
        for i=1:N                                                           %ypologizw to kostos thn lyshs poproekypse apo diaforopoihsh sto cost5
            for j=1:N
                cost5=cost5+standarC(i,j)*Xdiafor(i,j);   
            end
        end
        
        cost5
         
        if (cost5<cost)                                                     %kai an einai kalyterh apo thn hdh yparxousa thn kataxwrw ws nea beltisth
           X=Xdiafor;                                                           
           Xdiadromes=Xdiafordiadromes;
           cost=cost5;                                                     
        end
        
    end
    
    
end                                                                         %telos progrmmatos

cost

Xcordinates= Nxy(1,1);
Ycordinates= Nxy(1,2);
max=count;
startpoint=-1;
for i=1:count
    startpoint=1;
    acefound=1;
    while (acefound~=2)
        for j=1:N
            if (Xdiadromes(startpoint,j)==i)
          
                Xcordinates= [Xcordinates;Nxy(j,1)];
                Ycordinates= [Ycordinates;Nxy(j,2)];
                startpoint=j;
            
                if (startpoint==1)
                    if(i==1) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-r')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==2) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-g')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==3) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-b')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==4) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-y')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==5) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-m')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==6) && (i<max)
                        plot(Xcordinates,Ycordinates,'o-c')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i==7) && (i<=max)
                        plot(Xcordinates,Ycordinates,'o-k')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    if(i>=8)
                        plot(Xcordinates,Ycordinates,'o-r')
                        hold on
                        Xcordinates=Nxy(1,1);
                        Ycordinates=Nxy(1,2);
                        acefound=2;
                        break
                    end
                    
                end
            end
        end
                        if (i==max)
                            hold off
                            acefound=2;
                            break
                        end
    end
end


clientList=1;
start=1;
MAX=count;
for i=1:count
    acefound=1;
    start=1;
    while(acefound~=2)
        for j=1:N
            if(Xdiadromes(start,j)==i)
                clientList=[clientList;j];
                start=j;
                break
            end
        end
        if(start==1)
            clientList
            acefound=2;
            clientList=0;
            break
        end

    end
    
end