
clc;
clear;
close all;

%% Problem Definition

model=SelectModel();  % Select Model of the Problem
global c cycle_time;
cycle_time=1/3; %cycletme=20 minutes
c=12;
S=model.S;
s=1;
model.numeveihcle=[5 5 5 5];
model.Usedveihcle=model.numeveihcle;
it=1;
global directindex;
model.directindex=3.5;% max accaptable
initialindex=3;
directindex=initialindex;%%%in each
Pop=2;%initial pop for multi start SA
global next_helpstatue;
next_helpstatue=zeros(1,S);
%next_helpstatue=[0,1,0,0];
global before_helpstatue;
before_helpstatue=zeros(1,S);
%before_helpstatue=[0,0,1,0];
while s~=S+1

CostFunction2=@(q) MyCost2(q,model,s);
if (next_helpstatue (s)+before_helpstatue (s))~=0
    disp(['helping station' num2str(s)]);
else
    disp(['station number' num2str(s)]);
end
    %% SA Parameters
MaxIt=300;      % Maximum Number of Iterations                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
MaxIt2=250;      % Maximum Number of Inner Iterations 
T0=1;          % Initial Temperature
alpha=0.99;     % Temperature Damping Rate
%% Initialization
% Set Initial Temperature
T=T0;
% Array to Hold Best Cost Values
global BestCost;
BestCost(1:MaxIt,s)=zeros(MaxIt,1);
 
% Create Initial Solution
for p=1:Pop
   
   x(s,p).Position=CreateRandomSolution(model,model.Usedveihcle(1,s),s);
[x(s,p).Cost, x(s,p).Sol,x(s,p).Cost2 ]=CostFunction2(x(s,p).Position); 
end
% Update Best Solution Ever Found
BestSol(s)=x(s,1);
feasBestSol(s)=x(s,1);
%% SA Main Loop
for it=1:MaxIt
   sabet=model.directindex/initialindex;
      directindex=initialindex+((it^sabet))*(model.directindex-initialindex)/((MaxIt-100)^sabet);
   if (next_helpstatue(s)+before_helpstatue(s))>0 || it> MaxIt-100 
    directindex=model.directindex;   
   end
   if s==4
       directindex=2.7;
   end
      if it>300 && BestCost(it-1,s)==BestCost(it-300,s)
       break
      end
   [x(s,p).Cost, x(s,p).Sol,x(s,p).Cost2]=CostFunction2(x(s,p).Position);
    for it2=1:MaxIt2
        for p=1:Pop
        global good;
        good=1;
        
        % Create Neighbor
        while  good~=0  
        xnew.Position=CreateNeighbor(x(s,p).Position,model.I(s),model.Usedveihcle(s));
        [xnew.Cost, xnew.Sol,xnew.Cost2]=CostFunction2(xnew.Position);
        if (next_helpstatue(s)+before_helpstatue(s)==0)
            break
        end
        end
        
        if xnew.Cost<=x(s,p).Cost
            % xnew is better, so it is accepted
            x(s,p)=xnew;
            
        else
           
            % xnew is not better, so it is accepted conditionally
            delta=(xnew.Cost-x(s,p).Cost)/x(s,p).Cost;
            pi=exp(-delta/T);%-0.1*(1-x(s,p).Sol.IsFeasible);
            
            if rand<=pi
                x(s,p)=xnew;
            end
        end
            
        % Update Best Solution
        if x(s,p).Cost<=BestSol(s).Cost  %&&  x(s).Sol.IsFeasible==1
            BestSol(s)=x(s,p);
            if  x(s,p).Sol.IsFeasible==1
            feasBestSol(s)=x(s,p);
            
            end 
        end 
        end  
    end
    
    
    % Store Best Cost
    BestCost(it,s)=BestSol(s).Cost;
    
    % Display Iteration Information
    if BestSol(s).Sol.IsFeasible
        FLAG=' *';
    else
        FLAG='';
    end
    disp(['Iteration ' num2str(it) ':Best Cost= ' num2str(BestCost(it,s)) FLAG   '   -time: ' num2str(sum(BestSol(s).Sol.numetimefine)) '  -vehicles:' num2str(BestSol(s).Sol.Usedveihcle)]);%'  getoff'  num2str(BestSol(s).Sol.numgetoff)
    
    % Reduce Temperature
    T=alpha*T;    
   
end
    directindex=initialindex+it*(model.directindex-initialindex)/MaxIt;
    if  feasBestSol(s).Sol.IsFeasible
        if BestSol(s).Sol.IsFeasible
            PlotSolution2(BestSol(s).Sol,model,s)
        else
         PlotSolution2(feasBestSol(s).Sol,model,s) 
         BestSol(s)=feasBestSol(s);
        end
    else    
      %if model.numeveihcle(s)== BestSol(s).Sol.Usedveihcle
      if s>1 && BestSol(s-1).Sol.Usedveihcle < model.numeveihcle(s-1)
          model.Usedveihcle(s)=model.Usedveihcle(s)+1;
          model.numeveihcle(s-1)=model.numeveihcle(s-1)-1;
          before_helpstatue(s)=before_helpstatue(s)+1;
      else
          model.Usedveihcle(s)=model.Usedveihcle(s)+1;
          if s~=S
            model.numeveihcle(s+1)=model.numeveihcle(s+1)-1;
            next_helpstatue(s)=next_helpstatue(s)+1; 
          else
             disp(['***infeasible solution*** ']);
          end
      end
      %end
      % felan farzkardam istagah akhar az nakojabad mashin migireh

      s=s-1;
 
    end    
    s=s+1;
 
%   filename=['R station' num2str(s) 'timewindow' num2str(model.TimeWindow) 'ratio' num2str(10*model.directindex)];
%save(filename);
end
figure(2);
PlotSolution(model);
%%printing outouts to file
filename=['Results timewindow' num2str(model.TimeWindow) 'ratio' num2str(10*model.directindex)];
save(filename);

saveas(figure(2),['timewindow' num2str(model.TimeWindow), 'ratio' num2str(10*model.directindex),'passengers.jpg']);
saveas(figure(1),['timewindow' num2str(model.TimeWindow),'ratio' num2str(10*model.directindex),'gragh.jpg']);
%{
fileID = fopen('ResultsTW3ratio5.txt','wt');
formatSpec = 'number of stations: %1.0f\nspeed of vehicles(km/h): %1.2f\ncapacity of vehicles: %1.0f\nmaximum additional time ratio: %1.0f\n';
fprintf(fileID,formatSpec,model.S,model.v/1000,model.c,model.directindex);
for s=1:model.S
formatSpec = 'passangers station%1.0f: %1.0f';
fprintf(fileID,formatSpec,s,model.I(s));
end
for s=1:S
   formatSpec = '\nstation%1.0f X:';
   fprintf (fileID,formatSpec,s); 
   formatSpec = '%1.0f';
   fprintf (fileID,formatSpec,model.x0(1,s)); 
   formatSpec = '\nstation%1.0f Y:';
   fprintf (fileID,formatSpec,s); 
   formatSpec = '%1.0f';
   fprintf (fileID,formatSpec,model.y0(1,s)); 
end
for s=1:model.S
formatSpec = '\nBestSol station%1.0f arrival time:';
fprintf(fileID,formatSpec,s);
formatSpec = ' %1.4f';
fprintf(fileID,formatSpec,feasBestSol(s).Sol.AT);
end

fclose(fileID);
%}