%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  PDPTW for Multi-EV
%  Saman Ahmadi
%  15/08/2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
include "subcircuit.mzn";
include "nvalue_fn.mzn";

%%%%% Defining Required Sets of Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
string: Name; % Instance Name
int: Dist_Scale; % Distance Unit Scale
int: Time_Scale; % Time Unit Scale
float: Energy_Scale; % Energy Unit Scale
int: Chg_Scale; % Charger Points Duplicate Factor
int: Veh_Num; % Number of Vehicles
int: Veh_Type_Num; % Number of Vehicle Types
int: Chg_Type; % Number of Charger Types
array [1..Veh_Num] of int: Veh_Type; % Types of Vehicles in the Instance
array [1..Veh_Type_Num] of int: Veh_Capa; % Capacity of Each Vehicle Type
int: Veh_Speed; % Vehicle Speed
float: Inv_Veh_Speed; % Inverse Vehicle Speed
array [1..2] of int: Inst_Size; % Size of the 2-d Instance Matrix (Demand Details)
int: Dem_Num; % Number of {Pick-up, Delivery} Demand Pairs
int: Quick_Chg_Num; % Number of Quick Chargers
int: Fast_Chg_Num; % Number of Fast Chargers
int: Normal_Chg_Num; % Number of Normal Chargers
int: Chg_Num = Quick_Chg_Num + Fast_Chg_Num + Normal_Chg_Num; % Number of All Chargers
array [1..2] of float: Dechg_Rate;  % Average Decharge Rates alpha_3 and beta_3
array [1..Veh_Type_Num] of int: Energy_Min; % Minimum Energy for all Vehicle Types
array [1..Veh_Type_Num] of int: Energy_Max; % Maximum Energy for all Vehicle Types
array [1..Veh_Num] of int: Energy_Init; % Initial Energy Level of Vehicles in the Instance (in %)
int: Cur_Time; % Current/Start Time
int: Max_Time; % Upper Time Limit of the Instance
array [1..Veh_Type_Num,1..Chg_Type] of int: Max_Chg_Time; % Maxium Charging Time of Charger Types per Vehicle Type
int: Passenger_Weigth; % Passenger Weigth
float: Detour_Fac; % Maximum Detour Between Pickup and Dropoff
array [1..6,1..6] of float: Chg_Fcn; % Non-Linear Charger Function- Coefficints of the Linear Piecewise Charging Models
array [1..2,1..16] of float: Chg_Cut; % Non-Linear Charger Cutting Points, For time and Energy (for Charger Function)
array [1..Dem_Num,1..2] of int: Dem_Pair; % Demand Pair Matrix (pickup-dropoff) Locations
array [Locations,1..Inst_Size[2]] of float: Inst; % Details of Locations and Transport Requests in the instance
array [Locations,Locations] of int: Distance; % Distance Matrix
array [Locations,Locations] of float: alpha; % alpha Matrix (including speed profiles of links)
array [Locations,Locations] of float: beta; % beta Matrix (including speed profiles of links)
array [Locations,Locations] of float: alpha_avg;% alpha Average Matrix (using average coefficinets)
array [Locations,Locations] of float: beta_avg; % beta Average Matrix (using average coefficinets)

%%%%%%%%%%%%%%%%% Designing the Successor array and end points
% Succ=[pickup+dropoff+chargers+depot]
% order of chargers: Quick+Fast+Normal
% Setting endpoints in the Succ array
int: Dem_End=2*Dem_Num; % end of demand nodes: Number of Pickup + dropoff nodes
int: Chg_End=Dem_End+Chg_Num*Chg_Scale; % end of charger nodes
int: Veh_End=Dem_End+Chg_Num*Chg_Scale+Veh_Num; % End of depot nodes

%%%%%%%%%%%%%%%%% Creating useful sets
set of int: Vehicle = 1..Veh_Num; % All Vehicle Numbers
set of int: Locations = 1..Inst_Size[1]; % Set of All Locations in the Instance
set of int: All = 1..Veh_End; % Set of All Nodes (including copies of chargers)
set of int: Pickup = 1..Dem_Num; % Set of Pickup Nodes
set of int: Dropoff = Dem_Num+1..Dem_End; % Set of Dropoff Nodes
set of int: Dem_Range = 1..Dem_End; % Set of Demand (pickup+dropoff) Nodes
set of int: Dem_Chg_Range = 1..Chg_End; % Set of Demand+Charger Nodes (including copies of chargers)
set of int: All_Chg_Range = Dem_End+1..Chg_End; % Set of Charger Nodes (including copies of chargers)
set of int: Veh_Range = Chg_End+1..Veh_End; % Set of Vehicle Nodes
set of int: Quick_Chg_Range  = Dem_End+1..Dem_End+Quick_Chg_Num*Chg_Scale; % Set of Quick Charger Nodes 
set of int: Fast_Chg_Range = Quick_Chg_Range[length(Quick_Chg_Range)]+1..Quick_Chg_Range[length(Quick_Chg_Range)]+Fast_Chg_Num*Chg_Scale; % Set of Fast Charger Nodes
set of int: Normal_Chg_Range = Fast_Chg_Range[length(Fast_Chg_Range)]+1..Fast_Chg_Range[length(Fast_Chg_Range)]+Normal_Chg_Num*Chg_Scale; % Set of Normal Charger Nodes
set of int: Quick_Chg_Loc  =Dem_End+1..Dem_End+Quick_Chg_Num; % Set of Quick Charger locations in the Instance
set of int: Fast_Chg_Loc = Dem_End+Quick_Chg_Num+1..Dem_End+Fast_Chg_Num+Quick_Chg_Num; % Set of Fast Charger locations in the Instance
set of int: Normal_Chg_Loc = Dem_End+Fast_Chg_Num+Quick_Chg_Num+1..Dem_End+Fast_Chg_Num+Quick_Chg_Num+Normal_Chg_Num; % Set of Normal Charger Locations in the Instance

%%%%% Defining Required Constant Arrays %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% We now use the charger function and charger cutting points to creat the charging array per charger and vehicle type
%% The output array is Energy mapped based on the time unit
array [1..Veh_Type_Num,1..Chg_Type,0..max(Max_Chg_Time)*Time_Scale] of int: Energy_Chg=
array3d (1..Veh_Type_Num,1..Chg_Type,0..max(Max_Chg_Time)*Time_Scale,
[ round(if j<Chg_Cut[m,2+4*i]*Time_Scale then j*Chg_Fcn[3*(m-1)+i,1]*Energy_Scale/Time_Scale+Chg_Fcn[3*(m-1)+i,4]*Energy_Scale elseif j<Chg_Cut[m,3+4*i]*Time_Scale then j*Chg_Fcn[3*(m-1)+i,2]*Energy_Scale/Time_Scale+Chg_Fcn[3*(m-1)+i,5]*Energy_Scale elseif j<Chg_Cut[m,4+4*i]*Time_Scale then j*Chg_Fcn[3*(m-1)+i,3]*Energy_Scale/Time_Scale+Chg_Fcn[3*(m-1)+i,6]*Energy_Scale else Energy_Max[m]*Energy_Scale endif )| m in 1..Veh_Type_Num,i in 1..Chg_Type, j in 0..max(Max_Chg_Time)*Time_Scale] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Now we build energy arrays depending on the model
%%% ONLY ONE OF THE FOLLOWING MUST BE USED IN THE MODEL AT A TIME

%% For Base OR Mass only
% array [1..2,Locations,Locations,0..max(Veh_Capa)] of int: Dechg_Mat % Energy Matrix in the resolution, 1..Veh_Type_Num==1..1
%  = array4d(1..2,Locations,Locations,0..max(Veh_Capa),
%            [ myround(Energy_Scale*(Dechg_Rate[1]*m*Passenger_Weigth + Dechg_Rate[2])*Distance[i,j]) | t in 1..2, i,j in Locations, m in 0..max(Veh_Capa) ]);

%% For Gradient+Mass OR Gradient only
% array [1..2,Locations,Locations,0..max(Veh_Capa)] of int: Dechg_Mat % Energy Matrix in the resolution, 1..Veh_Type_Num==1..1
%  = array4d(1..2,Locations,Locations,0..max(Veh_Capa),
%            [ myround(Energy_Scale*(alpha_avg[i,j]*m*Passenger_Weigth + beta_avg[i,j])) | t in 1..2, i,j in Locations, m in 0..max(Veh_Capa) ]);

% For Gradient+Speed+Mass OR Gradient+Speed
array [1..2,Locations,Locations,0..max(Veh_Capa)] of int: Dechg_Mat % Energy Matrix in the resolution, 1..Veh_Type_Num==1..1
 = array4d(1..2,Locations,Locations,0..max(Veh_Capa),
           [ myround(Energy_Scale*(alpha[i,j]*m*Passenger_Weigth + beta[i,j])) | t in 1..2, i,j in Locations, m in 0..max(Veh_Capa) ]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Now we add mass into the model
%%% ONLY ONE OF THE FOLLOWING MUST BE IN THE MODEL AT A TIME

%% To Enable Mass in the model
constraint forall(i in All) (Link_Utl_Mass[i] = Link_Utl[i]);

%% To Disable Mass in the model
% constraint forall(i in All) (Link_Utl_ener[i] = 0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assistant constraints

%% To Disable Quick chargers
% constraint forall(i in Quick_Chg_Range)(Suc[i]=i);

%% To Disable Fast chargers
constraint forall(i in Fast_Chg_Range)(Suc[i]=i);

%% To Disable Normal chargers
constraint forall(i in Normal_Chg_Range)(Suc[i]=i);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Additional Arrays of constant

% Creating time array using speed and average speed
array [Locations,Locations] of int: Time= array2d(Locations,Locations,[round(Time_Scale*Inv_Veh_Speed*Distance[i,j]) | i,j in Locations]);

% Setting the initial energy in units of energy based on the initial % of energy and capacity
array [1..Veh_Num] of int: Energy_Start=[round(Energy_Init[i]/100*Energy_Scale*Energy_Max[Veh_Type[i]]) | i in 1..Veh_Num];

% Setting energy lower and upper bounds if set by user
array [1..Veh_Type_Num] of int: Energy_Low_Limit=[round(Energy_Min[i]/100*Energy_Scale*Energy_Max[i]) | i in 1..Veh_Type_Num];
array [1..Veh_Type_Num] of int: Energy_End_Limit=[round(100/100*Energy_Scale*Energy_Max[i]) | i in 1..Veh_Type_Num];

% Establishing the actual location of nodes
array [Vehicle] of int: Veh_Loc=[Inst_Size[1]| i in 1..Veh_Num]; % All vehicles at the same depot location
array [1..Chg_Num*Chg_Scale] of int: All_Chg_Loc=[i | i in Quick_Chg_Loc, j in 1..Chg_Scale]++[i | i in Fast_Chg_Loc, j in 1..Chg_Scale]++[i | i in Normal_Chg_Loc, j in 1..Chg_Scale]; % All of Charging Locations
array [Dem_Range] of int: Dem_Loc=[Dem_Pair[i,j] | j in 1..2, i in 1..Dem_Num]; % Demand Locations
array [All] of int: All_Loc=Dem_Loc++All_Chg_Loc++Veh_Loc; % All Actual Locations

% Establishing charger types of charger nodes
array [1..Chg_Num*Chg_Scale] of int: All_Chg_Type=[1 | i in Quick_Chg_Range]++[2 | i in Fast_Chg_Range]++[3 | i in Normal_Chg_Range];

% Establishing Delta Utilisation array for all nodes
array [All] of int: Del_Utl= [ myround(Inst[Dem_Loc[i],5]) | i in Dem_Range] ++ [ 0 | i in All_Chg_Range] ++ [ 0 | i in Veh_Range];

% Establishing Service Time array for all nodes
array [All] of int: Service_Time= [ round(Time_Scale*Inst[Dem_Loc[i],4]) | i in Dem_Range] ++ [ round(Time_Scale*Inst[All_Loc[i],4]) | i in All_Chg_Range] ++ [ 0 | i in Veh_Range] ;

%%%%% Time Window Configuration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
array [All] of int: Time_Arr_Low=[round(Time_Scale*Inst[All_Loc[j],6]) |j in All]; % Lower time window of nodes
array [All] of int: Time_Arr_Up=[round(Time_Scale*Inst[All_Loc[j],7]) |j in All]; % Max time of nodes

%%%%% Defining Required Variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
array [All] of var All: Suc; % Successor Representation in One Array
array [All] of var 0..Veh_Num: Veh_ID; % Vehicle ID is used in All nodes
array [All] of var Time_Scale*Cur_Time..Time_Scale*Max_Time: Time_Arr; % Arriving Time to All Nodes
array [All] of var 0..max(Veh_Capa): Link_Utl; % Utilisation in All nodes at departure
array [All] of var 0..max(Veh_Capa): Link_Utl_Mass; % Utilization in All nodes needed for Energy Matrix
array [All] of var 0..round(Energy_Max[2]*Energy_Scale): Energy_Arr; % Energy at Arrival for All Nodes
% array [1..Chg_Num*Chg_Scale] of var bool: Chg_ID; % Charger ID in Charger Stops Range
array [1..Chg_Num*Chg_Scale] of var 0..max(Max_Chg_Time)*Time_Scale: Chg_Time_Arr; % Arrival time offset for chargers
array [1..Chg_Num*Chg_Scale] of var 0..max(Max_Chg_Time)*Time_Scale: Chg_Time_Dep; % Departure time offset for chargers
array [1..Chg_Num*Chg_Scale] of var 0..max(Max_Chg_Time)*Time_Scale: Chg_Time; % Charging time for charger nodes

%% Assistant Variables
var int: Veh_Used=nvalue([Veh_ID[i] | i in Dem_Range]); % Number of Vehicle Used
var int: Quick_Chg_Used=sum(i in Quick_Chg_Range)(Suc[i]!=i); % Number of Quick Charger Used
var int: Fast_Chg_Used=sum(i in Fast_Chg_Range)(Suc[i]!=i); % Number of Fast Charger Used
var int: Normal_Chg_Used=sum(i in Normal_Chg_Range)(Suc[i]!=i); % Number of Normal Charger Used
var int: All_Chg_Used=sum(i in All_Chg_Range)(Suc[i]!=i); % Number of chargers visited
var int: Total_Time_Chg=sum(i in All_Chg_Range)(Chg_Time[i-Dem_End]); % Total Charging Time
var int: Total_Time_Drive=sum(i in All)(Time[All_Loc[i],All_Loc[Suc[i]]]); % Total Driving Time
var int: Total_Energy_Chg=sum(i in All_Chg_Range)(Energy_Arr[Suc[i]]-Energy_Arr[i]+Dechg_Mat[Veh_Type[Veh_ID[i]],All_Loc[i],All_Loc[Suc[i]],Link_Utl_Mass[i]]); % Calculating Total energy charged
var int: Total_Energy_Drive=sum(i in All) (Dechg_Mat[Veh_Type[Veh_ID[i]],All_Loc[i],All_Loc[Suc[i]],Link_Utl_Mass[i]]); % Total driving energy

%%%%% Adding Required Constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
constraint subcircuit(Suc); % Creating circuits/trips over nodes

%%%%% Applying Vehicle Constraints
constraint forall(i in 1..Veh_Num)(Veh_ID[Chg_End+i]=i); % Setting Vehicle ID for Vehicles in the depot
constraint forall(i in All where Suc[i]<=Chg_End)(Veh_ID[Suc[i]]=Veh_ID[i]); % Propagating Vehicle ID Inside Tours
constraint forall(i in 1..Dem_Num) (Veh_ID[i]=Veh_ID[Dem_Num+i]); % Pick-up and Delivery in the Same Vehicle
constraint forall(i in Veh_Range) (Suc[Dem_End+1]!=i); % No Depot Visit after pickup
constraint forall(i in 1..Veh_Num-1)(Suc[Veh_End-i]<Suc[Veh_End-i+1]); % Breaking symmetry for Vehicles at depot
constraint forall(i in Veh_Range)(Suc[Chg_End]!=i); % Breaking Symmetry: Last charger node is at the depot

%%%%% Applying Time Constraints
constraint forall(i in Veh_Range) (Time_Arr[Suc[i]]>=Cur_Time+Service_Time[i]+Time[All_Loc[i],All_Loc[Suc[i]]]); % Setting Departure Time for Vehicles in the Range
constraint forall(i in 1..Dem_Num) (Time_Arr[Dem_Num+i]>Time_Arr[i]); % Drop off Time After Pick-up Time
constraint forall(i in 1..Dem_Num) (Time_Arr[Dem_Num+i]-Time_Arr[i]<=Detour_Fac*Time[All_Loc[i],All_Loc[Dem_Num+i]]); % Limiting trip time for demands, Considering Detour Factor
constraint forall(i in Pickup)(Time_Arr[i]>=Time_Arr_Low[i]); % Respecting Lower Arrival Time at pickup nodes
constraint forall(i in Pickup)(Time_Arr[i]<=Time_Arr_Up[i]); % Respecting Upper Arrival Time at pickup nodes
constraint forall(i in Dem_Range)(Time_Arr[Suc[i]]>=Time_Arr[i]+Service_Time[i]+Time[All_Loc[i],All_Loc[Suc[i]]]); % Respecting Travel Time for non-charger nodes
constraint forall(i in All_Chg_Range where Suc[i]!=i)(Time_Arr[Suc[i]]>=Time_Arr[i]+Service_Time[i]+Chg_Time[i-Dem_End]+Time[All_Loc[i],All_Loc[Suc[i]]]); % Respecting Travel Time for charger nodes
constraint forall(i in All_Chg_Range)(Suc[i]==i <-> Time_Arr[i] = 0); % Fixing arrival time for unused charger nodes

%%%%% Applying Utility Constraints
constraint forall(i in Veh_Range) (Link_Utl[i]=0); % Setting Initial Utilization for Vehicles at depot
constraint forall(i in All_Chg_Range) (Link_Utl[i]=0); % No passenger while charging
constraint forall(i in All)(Link_Utl[Suc[i]]=Link_Utl[i]+Del_Utl[Suc[i]]); % Respecting vehicle utilisation at All Nodes
constraint forall(i in Dem_Range) (Suc[i]!=i); % All demands need to be served, no self loop for demand nodes

%%%%% Applying Energy Constraints
constraint forall(i in Veh_Range)(Energy_Arr[Suc[i]]=Energy_Start[Veh_ID[i]]-Dechg_Mat[Veh_Type[Veh_ID[i]],All_Loc[i],All_Loc[Suc[i]],Link_Utl_Mass[i]]); % Setting energy at departure
constraint forall(i in Dem_Range)(Energy_Arr[Suc[i]]=Energy_Arr[i]-Dechg_Mat[Veh_Type[Veh_ID[i]],All_Loc[i],All_Loc[Suc[i]],Link_Utl_Mass[i]]); % Respecting energy at non-charger nodes
constraint forall(i in All_Chg_Range)((Energy_Arr[Suc[i]])=Energy_Chg[Veh_Type[Veh_ID[i]],All_Chg_Type[i-Dem_End],Chg_Time_Dep[i-Dem_End]]-Dechg_Mat[Veh_Type[Veh_ID[i]],All_Loc[i],All_Loc[Suc[i]],Link_Utl_Mass[i]]); % Respecting energy at charger nodes

%%%%% Applying Charging Constraints
constraint forall(i in Pickup) (Suc[i]<=Dem_End); % No charge after pickup
constraint forall(i in All_Chg_Range) (Suc[i]<=Dem_Num \/ Suc[i]>Chg_End \/ Suc[i]=i); % No dropoff/recharge after charge
constraint forall(i in Veh_Range)(Suc[i]<=Dem_End); % No charging visit after leaving the depot
constraint (Total_Energy_Chg>=0); % Prevent negative energy at chager nodes
constraint (All_Chg_Used <= Veh_Num); % Only one charging visit per vehicle
constraint forall(i,j in All_Chg_Range where i<j) ( (Suc[i]!=i /\ Suc[j]!=j) -> ( Veh_ID[i]!=Veh_ID[j] )); % Only one charging visit per vehicle
% constraint forall(i in All_Chg_Range)(Suc[i]=i <-> Chg_Time_Dep[i-Dem_End] = 0); % Unused charger, fix charging time offset
constraint forall(i in All_Chg_Range)(Chg_Time[i-Dem_End]=Chg_Time_Dep[i-Dem_End]-Chg_Time_Arr[i-Dem_End]); % Calculating charging time per charger node
constraint forall(i in All_Chg_Range)(Chg_Time_Dep[i-Dem_End]>=Chg_Time_Arr[i-Dem_End]); % Prevent negative charging times

%%%%% Binding/Tuning the start charge time offset to the closest integer
constraint forall(i in All_Chg_Range)(Energy_Arr[i]>=Energy_Chg[Veh_Type[Veh_ID[i]],All_Chg_Type[i-Dem_End],Chg_Time_Arr[i-Dem_End]]); 
constraint forall(i in All_Chg_Range)(Energy_Arr[i]<Energy_Chg[Veh_Type[Veh_ID[i]],All_Chg_Type[i-Dem_End],Chg_Time_Arr[i-Dem_End]+1]); 

%%%%% Set upper bounds on charging times depending on the charger type and max charging time
constraint forall(i in All_Chg_Range)(Chg_Time_Arr[i-Dem_End] <= Time_Scale*Max_Chg_Time[Veh_Type[Veh_ID[i]],All_Chg_Type[i-Dem_End]]);
constraint forall(i in All_Chg_Range)(Chg_Time_Dep[i-Dem_End] <= Time_Scale*Max_Chg_Time[Veh_Type[Veh_ID[i]],All_Chg_Type[i-Dem_End]]);
constraint forall(i in All_Chg_Range)(Chg_Time[i-Dem_End] <= Time_Scale*Max_Chg_Time[Veh_Type[Veh_ID[i]],All_Chg_Type[i-Dem_End]]);

%%%%% Objective Function and Search %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
solve 
:: int_search(Suc, input_order, indomain_min, complete)
  minimize ((Total_Time_Drive+Total_Time_Chg)); % Mininizing Total travel time and charging time

%%%%% Output and Display %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
output [  "Instance Name      = ",show(Name)++ "\n"++
          "Total Drive Time   = ",show(Total_Time_Drive/Time_Scale)++ " min\n"++
          "Total Charge Time  = ",show(Total_Time_Chg/Time_Scale)++ " min\n"++                 
          "Tot Driving Energy = ",show(Total_Energy_Drive/Energy_Scale/1000)++ " kWh\n"++
          "Tot Charged Energy = ",show(Total_Energy_Chg/Energy_Scale/1000)++ " kWh\n"++
          "Quick Charger Used = ",show(Quick_Chg_Used)++ "\n"++
          "Fast Charger Used  = ",show(Fast_Chg_Used)
%           "Location Array    = ",show([All_Loc[Suc[i]]|i in All])
];
%%%%% Additional Functions Definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function int: myround(float: c) =
if c>=0 then round(c) else -round(-c) endif;