A HETEROGENEOUS ELECTRIC TAXI FLEET ROUTING PROBLEM WITH RECHARGING STATIONS TO MAXIMIZE THE COMPANY’S PROFIT

. During the past years, many kinds of research have been done in order to reduce the cost of transportation by using different models of the vehicle routing problem. The increase in the amount of pollution caused by vehicles and environmental concerns about the emission of greenhouse gases has led to the use of green vehicles such as electric vehicles in the urban transport fleet. The main challenge in using electric vehicles with limited battery capacity is their long recharging time. For this purpose, several recharging stations are considered in the transportation network so that if the battery needs to be recharged, the electric vehicle can recharge and complete its journey. On the other hand, due to the limited amount of the electric vehicle’s energy, the fuel consumption of this fleet is highly dependent on their load, and it is necessary to consider their load in the planning. In this article, the problem of routing electric taxis is presented considering the economic and environmental aspects of implementing electric taxis for city services. Despite other studies that have only focused on reducing energy consumption or minimizing distance traveled by electric vehicles, for the first time, the problem of urban electric taxi routing has been modeled by considering different types of electric taxis with the aim of achieving the maximum profit of this business. The use of a heterogeneous fleet in this study leads to wider coverage of different types of demand. Therefore, a mathematical programming model is presented to formulate the problem. Then, several problem examples are designed and solved for validation purposes, and the simulated annealing algorithm (SA) will be introduced and used to solve large-scale problems.


Introduction
According to the United Nations, economic growth has increased the demand for transportation, with 44% and 78% of goods and passengers being transported by land, respectively, which has attracted much attention to the transportation system.With the evolution of technology, logistic and supply chain management practices have developed and subjects such as greening, globalization, agility and automation [23].In the field of logistics, green logistics has attracted more attention as an example of sustainable development.The sustainable perspective has encouraged governments and industries to enhance environmental sustainability and become more ecofriendly by improving their activities [29,35].Transportation is the biggest source of pollution in logistics.In Canada, for example, transportation accounted for 27% of greenhouse gas emissions in 2007, and in the United States, transportation accounted for 28% of national greenhouse gas emissions [38].
Fossil energy is the largest source of carbon dioxide and greenhouse emissions, and one of the main global issues is decreasing greenhouse gas emissions [6].The U.S. Environmental Protection Agency persuades suppliers and buyers to produce green products to decrease carbon emission [1].Research shows that people are drawn to non-fossil fuels.Environmental protection and optimal use of non-renewable resources are the priority of societies in today's world [20].To achieve these goals, it is necessary to minimize the use of fossil fuels in transportation.Experts suggest different methods to achieve this goal.For example, minimizing the distance traveled by using automated vehicles that travel the shortest route and have scheduled stops reduces the use of fossil fuels [10].The introduction of green vehicles, such as hydrogen vehicles and electric vehicles, is another good example of communities struggling to use non-fossil fuels [8].The environmental benefits of electric vehicles, including reduced transportation costs and environmental detrimental effects, have increased the logistics community's interest in incorporating electric vehicles into its transport fleet.However, limitations such as long charging times and limited access to charging facilities make using electric vehicles challenging.Due to technological limitations, electric vehicles usually have a short range of 100-150 miles.They need refueling or recharging along the way.The number of recharge stations is very limited.For example, there are only 1626 recharge stations across Canada [13].
The problem of electric vehicle routing is in fact, a development of the traditional vehicle routing problem, which deals specifically with finding optimal routes for electric vehicles, taking into account battery limitations and charging operations [32].To obtain more realistic results for the problem of electric vehicle navigation, the calculation of energy consumption must take into account several aspects such as road conditions, vehicle technical characteristics, vehicle load and environmental conditions.But the main challenge is considering any additional energy consumption aspects that complicate the calculations.Therefore, the various aspects considered for calculating energy consumption affect the acceptable calculation time and the degree to which the solution can be implemented in practice [16].
In existing studies, energy consumption calculations can be divided into linear deterministic functions, nonlinear deterministic functions and stochastic functions.Linear deterministic functions can be classified based on the parameters of the equation used into five subgroups: distance traveled, vehicle speed, vehicle load, road slope, and time.Linear equations have been used more by researchers to formulate the problem of electric vehicle routing as a complex integer programming model.If the equation is based on distance, energy consumption is determined by the energy consumption scale for a given distance.Similarly, other linear functions use a set of parameters to calculate energy consumption [24].
In this study, for the first time, the routing problem of the urban electric taxi company has been investigated in terms of the profitability of this business and the economic justification of its implementation.Considering the environmental benefits of its implementation, the establishment of different types of electric taxis in the urban transport fleet has an economic justification if it is profitable.Driven by the characteristics of electric vehicles mentioned before, efficient routing management of these vehicles is essential.The company has a charging station in the electric taxi depot with a limited number of charging stations in the city.If the remaining charge of the electric taxi is not enough to make a trip, it can perform charging operations during the trip by visiting the charging stations.In our article, we describe this problem with a mathematical model, which incorporates the interests of the electric taxi company in the objective function.
The rest of this paper is organized as follows.In Section 2, a literature review is presented.In Section 3, the problem is described, and a mathematical model is provided.In Section 4, our solution algorithm is introduced and explained.Section 5, describes the experiments and model validation.Section 6, shows the results of numerical experiments and insights.Finally, Section 7 concludes our study and mentions important results and future opportunities.

Related works
As the main issue in the planning and operation of service systems, procurement, and distribution of goods, the vehicle routing problem is of great theoretical importance and scientific value, and there are many studies in this field.In this regard, the routing of electric vehicles as a newer version of fossil vehicles has received more attention in recent years.
Paz et al. [25] studied both conductive charging and battery replacement methods for electric vehicles.The authors introduced three different mathematical formulas for three strategies, one for conductive charging, in which electric vehicles are charged at customer locations or conventional charging stations, and one for battery replacement, where electric vehicles can be charged at charging stations.Normally charged, and finally, for both conductive charging and battery replacement options where conventional charging stations are used to replace batteries while customers location are used for conductive charging.Meng and Ma [21] presented an optimization problem for routing electric vehicles with time windows based on two charging methods and solved their model using the ant colony algorithm.Sayarshad et al. [31] studied the problem of dynamic routing of electric taxis with battery replacement stations with a predictive policy using the Markov decision-making process to allocate the fleet of electric taxis to customers on the assumption of elastic demand.Jie et al. [12] introduced the problem of electric vehicle routing as a two-stage optimization problem and proposed a new algorithm for optimizing ant colonies at two-stage capacity vehicle routing and the charge of vehicles for fixed-route travel.Zhou and Tan [41] simultaneously examined the location of battery replacement stations and the routing of electric vehicles.In their study, battery replacement station locations and electrical vehicle routing problems are optimized for scheduling material handling operations at vehicle assembly plants.
Nonlinear definite functions can be studied in heuristic approaches to more realistically simulate energy consumption.In this regard, Goeke and Schneider [9] introduced a large adaptive large neighborhood search algorithm that uses nonlinear formulas to determine energy consumption values.The authors tested the performance of this algorithm on a benchmark data set in which the problem size varied from 10 to 100 client nodes.Similarly, Preis et al. [27] proposed an adapted tabu search algorithm that considers nonlinear equations.The results of their studies show that nonlinear equations can be effectively applied to the problem of real electric vehicle routing using a heuristic approach.In more recent studies, nonlinear deterministic functions have been introduced in the variable neighborhood search algorithm.Zhang et al. [39] proposed an ant colony optimization approach considering nonlinear target functions.The authors compared the performance of their algorithm with the large adaptive neighborhood search algorithm proposed by Goeke and Schneider [9] and showed that the ant colony optimization algorithm could provide better results for the routing of electric vehicles in computational times.Li et al. [17] introduced a mixed-integer mathematical programming model considering the heterogeneous electrical fleet with time windows and simultaneous collection and delivery, aiming to minimize distance.In their study, energy consumption is proportional to the travel time, and a load of electric vehicles and the variable neighborhood search algorithm are used to solve the model.Li et al. [18] developed a comprehensive model for the electric vehicle routing problem, which uses speed, load and distance indices to calculate energy consumption and carbon emissions, then they create adaptive genetic algorithms based on optimizing hill climbing and neighborhood search.A neighborhood search strategy is used to obtain the optimal solution to satisfy battery life limitations and battery replacement stations.Lu and Wang [19] proposed an optimization algorithm based on two strategies along with dynamic capacity consideration.Finally, according to the characteristics of the electric vehicle routing problem, they have designed two search strategies for display and scheduling.Kim et al. [15] presented a robust electric vehicle routing model under road traffic flow uncertainty and used a genetic algorithm with three-layer chromosome encoding schemes to solve the proposed model.
In addition to definite functions, stochastic approaches have been used to estimate the energy consumption of electric vehicles.Pelletier et al. [26] presented the electric vehicle routing problem with a robust optimization framework to consider energy uncertainty and developed a two-step heuristic method based on large neighborhood search to solve large-scale models.Reyes-Rubiano et al. [28] investigated the Green Electric Vehicle Routing Problem.Their study considers stochastic travel times and driving range limitations when running the car battery.In order to design reliable Routing plans, a sim heuristic algorithm that combines Monte Carlo simulation with a multi-start metaheuristic is proposed.Their main goal is to create reliable, low-risk solutions to request remedial action.They also represent how using appropriate energy safety storage levels in the routing design phase can increase the reliability of distribution plans and reduce the total expected costs.Soysal et al. [34] developed a nonlinear programming model of mixed integers, assuming that the batteries of electric vehicles are randomly discharged.Considering uncertainty makes it possible to meet the desired level of service and helps reduce driver anxiety.Zhang et al. [40] studied the electric vehicles routing problem in a fuzzy environment by considering time windows, the possibility of partial charging, and uncertainty of service times, battery consumption, and travel times.A large adaptive neighborhood search algorithm with a fuzzy simulation method has been developed to solve the proposed model.Zhou and Zhao [42] proposed a hybrid electric vehicle routing problem considering battery swap and mixed time window constraints.To solve this problem, the authors developed a multi-objective whale optimization algorithm enhanced by a particle filters and Levy flights.Vahedi-Nouri et al. [36] investigated a shared capacity electric vehicle routing problem with the goal of cooperation to reduce operating costs and improve customer service levels.In this regard, two new bi-objective mathematical models have been developed for the problem under cooperative and non-cooperative strategies.
According to existing studies, the following research gaps exist in the literature: -In almost all studies, the amount of energy consumed by an electric vehicle for a given distance is determined using constant energy consumption per unit distance.However, in real driving conditions, the energy consumption of an electric vehicle is affected by various conditions, such as road conditions, traffic, weather, driver performance, and vehicle load.In addition, the energy consumption of the heating or air conditioning systems of electric vehicles can be analyzed according to the weather conditions.Because optimization algorithms are ultimately optimized, some routes will have a near-zero battery level when they reach a charging station or warehouse.Therefore, ignoring or simplifying more than the energy consumption function may lead to practically impossible paths to implement.-Most electric vehicle navigation research minimizes distance-based performance metrics, while few consider the total energy consumed or the number of charging stations used.The need to address the concepts of green transportation has led to the fact that the objective functions under consideration are more relevant to energy consumption or charging operations.However, the real goal in most applications is to minimize the total cost, which includes energy consumption, driver cost, and vehicle cost.-Heterogeneous fleets increase the complexity of the problem for both precision solvers and exploratory algorithms.Choosing the type of vehicle for a route directly affects the cost and time of travel or may lead to violations of constraints such as capacity or other related constraints [11].On the other hand, a heterogeneous fleet has the potential to reduce the total cost of transportation or energy consumption by selecting more suitable vehicles for the routes but has been less studied.
Our objective is to specifically support electric vehicles' use for urban operations.The use of electric taxis in the urban transport fleet can be implemented if it is profitable for governments and non-governmental companies.In this study, the amount of energy used is calculated by considering factors such as speed, load, road friction and distance traveled.Also, considering the types of electric taxis (economical, luxury, and vans), leads to answering a wide range of demands.In summary, the main contributions of this study are to model and solve a practical problem resulting from the fact that electric taxis must be able to perform their full travel routes in urban areas with their departing battery charge.Therefore, this study focuses on implementing the electric vehicle design in the urban fleet with the goal of profitability, despite other studies focusing only on reducing energy consumption or minimizing the distance traveled by electric vehicles.

Problem description and assumptions
In this section, we describe the mathematical model of electric taxis routing.Each route services a set of customer nodes and starts and ends at a given depot node.The problem aims to find the best route plan for electric vehicles that maximize company profits while satisfying operational restrictions for electric taxis.The basic assumptions of this study are summarized below: -Each trip begins and ends at the Depot node.
-Each customer node must be serviced exactly by an electric taxi.
-Each charging station can be visited by more than one electric vehicle.
-The location of the charging stations already exists and the distance traveled from each node to each charging station are specified.-When visiting the charging station and depot, the battery of the electric taxi is fully charged.
-Three types of electric taxis (economy, luxury and van) are provided for more service flexibility.
-Mixed integer model including battery limits, load capacity developed.
The calculation of battery energy consumption is required to formulate the electric vehicle routing model.The forces acting on the vehicle are due to internal resistance, tires and air.The total tensile force,   , can be estimated by the following equation: is the frontal surface area of the electric taxi (m 2 ),  is vehicle's speed (m/s),  mass of vehicle (Kg),  is the gravitational constant whose value is considered to be 9.81(m/v 2 ),  is the air density (kg/m 3 ),   and   show the rolling resistance coefficient and the drag resistance coefficients, which are determined from the experiment.A typical value for   is 0.015 and   for cars vary, a value of 0.3 is commonly used.The mechanical power (MP) of electric taxis is required to calculate energy consumption.MP can be determined from the drag force given above, road inclination angle, and vehicle velocity [2].
where  is the road inclination angle, and  represents acceleration (m/v 2 ).The calculation of the mechanical power required in each arc (, ) for an electric taxi of type  is as follows [4]: We assume that the parameters remain constant in each arc but may be different in the other arc.The electric taxi travels on an arc (, ) at an average speed of    and carries a load of   +    where    is the total cargo weight and   is the curb weight of the electric taxi .Traveling distance is indicated as   .  is the specific coefficient of the arc and  is the specific coefficient of the vehicle and are calculated as follows: Considering an electric taxi's motor efficiency, the mechanical energy required in each arc is converted to the electrical energy required by the car's electric motor [39].The battery energy required to move in the arc (, ) is calculated as follows: where   = 1.25 is motor efficiency and   = 1.11 is battery discharge efficiency [24].The carbon dioxide emissions of an electric taxi in each arc are estimated as follows, taking into account the above formula, battery charging efficiency, and rate of conversion of electrical energy into greenhouse gases.
Assume that the fuel required by the power plant is obtained from burning coal, then the rate of conversion of electrical energy into greenhouse gases is   = 0.69 (kg/kwh) and the recharging efficiency of a battery is   = 1.25 [8].For simplicity, we define ET as follows:

Notations and model formulation
In this subsection, a mathematical model of this problem is proposed.The notations used in the model are defined in Table 1.The mathematical formulation of the electric taxis routing problem can be stated as follows.
Mathematical formulation of the electric taxis routing problem can be stated as follows [32,33,39]: Subject to ∑︁ The objective function (9) aims to maximize the total net cash flow of the company in time periods by calculating the company's revenue and reducing energy consumption costs.The first part of the objective function shows the company's revenue from providing taxi services, which is calculated by deducting the amount received from customers' and taxi drivers' salaries.The second part of the objective function calculates the cost of energy consumption according to the formula provided for calculating energy consumption.
Constraint (10) ensures that each customer is visited only once.Constraint (11) ensure that each charging station and its associated dummy vertices have a maximum of one surrogate, and taxis can travel from each charging station to a maximum of one more vertex.Constraints (12) ensure that the number of input and output arcs to the customer nodes and charging stations is equal, ensuring the path's continuity.Constraint (13) ensures that any electric vehicle can only be used in one route plan.Constraints ( 14) and ( 15) adjusts the battery level after visiting customers and charging stations, also ensures that the battery level does not fall below the minimum acceptable value and constraints (16).Controls the battery level when moving from Depot nodes.Constraint (17) ensure that the weight is reduced after the customer arrives.Constraint (18) limits the maximum load carried by an electric taxi k and force if    = 0 then    = 0. Constraints (19) and ( 20) cover arrival and end times at any time after customers leave and cover it after leaving the charging station, depending on the time to charge the battery.Arrival times within time windows are guaranteed by constraints (21) for both customers and charging station nodes.Constraints ( 22) and ( 23) define the acceptable values for    and    .

Solution approach
Since the issue of vehicle routing is one of the Np-hard issues, so the model of this research is also part of this category.Despite the efficiency of accurate algorithms in solving difficult small-scale optimization problems, solving large-scale problems with these algorithms is very time-consuming and tedious [30].Unlike accurate methods, metaheuristic methods produce appropriate solutions in an acceptable time for large-scale problems.As the solution method to the electric vehicle routing problem, a simulated annealing (SA) algorithm is proposed.

SA algorithm-based meta-heuristic
SA is a local search-based algorithm for searching the answer space.This method is obtained from the thermodynamic process of cooling molten metals to obtain the lowest energy state [22].The main advantage of SA over random search methods is its ability to avoid falling into local optimizations when searching for local minimums.In the SA search strategy, the algorithm starts from an initial answer.At each stage, new answers are generated by random gestures in the neighborhood of the current answer.If the new answer improves the value of the objective function, it replaces the current answer, otherwise, the new answer has a certain chance of being accepted.A new answer is accepted if this random number is less than probability.At the beginning of the algorithm, worse answers are more likely to be accepted.This allows us to search for a larger area of the answer space.The temperature in each step remains constant until a constant number of repetitions is performed.It is then gradually reduced so that only effective moves are accepted in the end.The algorithm continues until the final temperature is reached [7].

The answer display structure of the SA algorithm
To display the answer, assume that the number of customers is equal to , the number of fuel stations is equal to FS, and the number of vehicles is equal to .The answer chromosome is a permutation with dimensions of 1 × ( + FS +  − 1).For example, suppose we have 2 electric taxis, 10 customers, and 3 charging stations, with a virtual node for each charging station, which means we have 6 charging stations and 10 customers.The answer chromosome is a 1 × 17 matrix (Fig. 1).On the answer chromosome, the numbers 1 to 10 indicate customers, and the numbers 11 to 16 indicate charging stations.It is also the number 17 separator.

Creating a neighborhood
In the refrigeration simulation algorithm, in order to reach a new solution, it is necessary to change the structure of the current solution to find a new neighborhood in the solution space.In this algorithm, we apply 3 neighborhood creation operators that are randomly selected in each iteration of the refrigeration simulation algorithm [37].4.1.2.1.Swap.In this operator, first 2 nodes are randomly selected and their places are changed with each other.In Figure 3, the Swap operator is applied to nodes 4 and 7. 4.1.2.2.Reversion.In this operator, two nodes are randomly selected and in addition to replacement, the nodes between them are reversed.The mechanism of this operator based on the assumption of selecting nodes 4 and 7 on chromosomes in Figure 4 is as follows.4.1.2.3.Insertion.In this operator, two nodes are randomly selected and the first node is removed and moved to the second node.The mechanism of Insertion operator based on the assumption of selecting nodes 4 and 7 on chromosomes in Figure 5 is as follows.

Numerical examples
In the following, first we run the proposed model with GAMS 24.1.2software to validate it.The proposed SA algorithm is coded in MATLAB programming language and is implemented on a computer with CPU Intel Core i7 and 8 GB of RAM.The optimal solution for small and medium-sized problems is obtained.Then we run the problems in larger sizes with 9, 11 and 13 nodes using the SA algorithm and examine the deviation of the algorithm's answer from the exact solution.

Computational experiments
As seen in the model presented in Section 3, this study can be applied to a wide range of network topologies in terms of the location of customers, charging stations and any battery capacity.Clearly, the relative position between customers, charging stations, and battery capacity strongly influence the routing strategies of electric taxis.
All small solved problems have only one depot located in the center of a 100 × 100 km grid.Customers and charging stations are randomly scattered in the network.Arc-specific coefficient (  ) is considered constant and is equal to 0.0981.Between any two locations in the grid, an electric taxi moves at a constant speed randomly selected from the uniform (20,80).We assume three types of electric taxis (economy, luxury and van) with curb weights of 8.5 kg, 9 kg and 15 kg.Table 2 shows the rest of the parameters.Travel time from node  to node  (  ) is logically derived from the relationship between speed and distance (  =   /  ).
It should be noted that different methods due to their nature prepare different parameters.example, battery capacity and average speed were calculated from catalogs, experts' opinions and their performance in past periods.Data such as income and expenses are obtained through communication with management and consultation with them.In addition, parameters such as service time in nodes and travel times are obtained through timing methods such as stopwatches and so on.
As the number of customers and charging stations increases, the dimensions of the problem increase.Solving large-scale problems by exact solution methods is not justified due to the exponential increase in solution time, so we solve them with the SA algorithm.

Parameters tuning
In this section, the Parameters selected for adjustment include external iteration count (MaxIt), internal iteration count (MaxSubIt), initial temperature (T0), and cooling rate (Alpha).These parameters were examined at three levels: low, medium and high, and the values of each parameter at each level were considered using other research.Table 3 presents the proposed levels for the parameters of the SA algorithm.
After calculating the design of experiments by the Taguchi method in MINITAB software, the optimal value of the parameters is as shown in Table 4. Taguchi method and experiment designs are statistical methods by which it is possible to create experimental designs that are not affected by sample selection or process.The Taguchi method is the most powerful method available to reduce product costs, improve quality and at the same time reduce the development gap.By implementing the L9 test plan on a sample problem and de-scaling the obtained objective function values and in order to minimize the de-scaled values, the signal-to-noise ratio diagram is obtained according to Figure 6.

Model validation
Validation investigates whether the outputs of the presented model are quantitatively and qualitatively consistent with independent observations of the real world and whether its processes have a range of correctness and consistency from the point of view of the model's application.In order to guarantee the validity and reliability of the model presented in Section 3, we examine the effect of payload capacity of taxis and the demand on the model.
In the following, two examples are presented with the aim of analyzing the impact of reducing payload capacity and increasing demand separately on the objective function and optimal routes.Tables 5-7 show the results of the changes of these two parameters on the model.In Table 5, a problem example with 6 nodes (3 customers, 2 charging stations and 1 depot node) is solved.
In Table 6, the payload capacity of taxis decreases, but the amount of customer demand remains constant.It is clear that the objective function is not improved by reducing the payload capacity and some optimal routes are changed, but total cargo weights remain constant.Also, due to the reduction of the payload capacity, the optimal route of electric taxi type 2 changes and needs to visit charging station 5. Now it is assumed that the payload capacity of taxis will remain constant and the amount of customer demand will increase.In Table 7, the amount of customer demand increases; consequently, the total cargo weights and the objective function decrease.As a result of increasing the total cargo weights, electric taxi type 1 covers the customer demand of type 3, and the optimal routes of taxis do not change.Hence, the results demonstrate an expected trend, which shows the validation of the proposed model.

Results and discussions
this section, the results of solving sample problems using the CPLEX solver in GAMS and the SA algorithm in MATLAB are presented.

Results
In Table 8, the results of solving sample problems by GAMS and SA algorithms are given.As it is evident from the results of Table 8, the SA algorithm spends less time than the exact solutions.Also, in small dimensions, the optimality gap between the SA algorithm and GAMS is small, indicating the SA algorithm's proper efficiency.
In Figure 7, it is clear that with the increase in the dimensions of the problem, the solution time in the exact method increases exponentially.According to Table 8, the gap between the exact solution method and meta-heuristic is below 4%, so the SA algorithm has a good performance and can be used to solve large-sized problems.
Table 9 shows the results of solving another sample problem and the optimal routes of electric taxis in each problem by the SA algorithm.As it is evident from the results of Table 8, it reaches the optimal and acceptable solution at the right time.Table 10 show the convergence of the SA algorithm.In problem 1, it is clear that the algorithm has reached convergence from the initial iterations.In problem 2, the algorithm has reached convergence from iteration 30 and in problem 3 from iteration 250 onwards.
In Table 11, we solved problems with larger dimensions and more flexibility with the SA algorithm.In these problems, unlike Table 9, where the issues were limited to 3 types of electric taxis and 2 charging stations, the number of electric taxis changes according to the number of charging stations.
Figures 8 and 9 show the graphical representation of problems 4 and 5.The coloring of the lines indicates the optimal routes of all types of electric taxis.

Sensitive analysis
Sensitivity analysis is used to apply changes in model parameters and their effect on the value of the objective function.For this purpose, we consider the battery capacity parameter of electric taxis, the demand and payload capacity of electric taxis.The range of parameter changes is from +20% to −20% with a step of 10%.
Figure 10 shows the changes in customer demand on the value of the objective function.With the increase in demand, the amount of load transported between the nodes will increase, increasing energy consumption and as expected, the increase in demand leads to a decrease in profits.Figure 11 shows the sensitivity analysis of the payload capacity of electric taxis.Due to the fact that the fleet of electric taxis is heterogeneous, reducing their capacity leads to a decrease in the profit of the electric taxis company.Therefore, it is recommended to the management to increase the payload capacity to achieve more profit.Also, according to Figure 10, the amount of profit decreases with the increase in demand, and the management can compensate for this decrease in profit by increasing the payload capacity.
As shown in Figure 12, the objective function increases with the increase of the battery capacity of electric taxis.Increasing the battery capacity may lead to fewer charging stations being visited, in which case the distance traveled will decrease and lead to a reduction in costs.Therefore, the management can also increase the battery capacity of electric taxis to achieve more profit.
Figure 13 shows the sensitivity of the objective function to changes in the unit cost of energy consumption, that the objective function decrease with the increase of this parameter.Also, Figure 14 shows the sensitivity analysis of all the parameters examined in this section on the objective function in a graph.

Discussion and managerial insights
In this study, the electric taxi routing problem has been modeled by considering the types of electric taxis (economy, luxury, and van) to meet a wider range of demands.Several limitations related to vehicle capacity, battery capacity, recharging operations, and energy consumption are considered in this paper.Due to the size of their batteries occupying a large amount of available space, electric vehicles have a more limited cargo volume than conventional vehicles; therefore, to ensure that justified answers are obtained in the routing system, this constraint was considered in the model.To bring the issue closer to the real world, this study focuses on the implementation of the electric vehicle design in the urban fleet, which aims to make the use of the electric taxi fleet profitable.In most studies, the main goal is to reduce fuel consumption or minimize the total distance traveled by the electric vehicle, and few studies have addressed the economic aspects of implementing an electric taxi fleet [3].Most studies consider a homogeneous fleet, and few studies have considered different types of electric vehicles with different characteristics to provide better services to different types of demand [12].In this study, the amount of fuel consumption depends on the amount of load carried by each type of taxi, the weight of different types of taxis, and quantities such as acceleration and speed that make the calculation of energy consumption closer to reality.This type of energy consumption calculation has been used in fewer studies to avoid nonlinear functions that lead to model complexity [5,41].Using a fleet of electric taxis for intra-city transportation, in addition to reducing greenhouse gas emissions, also has an economic justification and is profitable for private and public companies.To achieve more profit, it is recommended to increase the payload capacity.If the amount of demand increases and leads to a decrease in profit, the management can compensate for this decrease in profit by increasing the payload capacity.Increasing the battery capacity of electric taxis also makes more profit.

Conclusions
In this study, for the first time, a vehicle routing problem for urban electric taxis has been modeled considering the heterogeneous fleet and with the aim of achieving the maximum profit for this business, and the presented model was solved using the exact method and meta-heuristic SA algorithm, and the efficiency of the SA algorithm was shown in large-scale problems.Then, several sensitivity analyses were performed on the model's key parameters.
In future studies, considering a mixed fleet including electric and hybrid vehicles, applying other restrictions related to vehicle volume, such as determining the arrangement of goods and limiting the tolerable pressure of items, are some of the innovations that can be developed on this model.Although this study is modeled for an intra-city transport fleet, it can also be extended to out-of-town trips.To get closer to reality, the partial charging policy can be used to charge the batteries of electric taxis.Mobile charging stations can be considered in the models, thus making transportation more flexible.Another research suggestion is to consider the possibility of battery replacement in charging stations and the cost balance between battery replacement and charging.Other innovative algorithms can also be used to determine charging station visits and generate an initial quality response to improve the proposed solution algorithm.Finally, with the increasing development of driverless vehicles and artificial intelligence, more problems, models and algorithms can be investigated in relation to the routing of electric vehicles.

Figure 1 .
Figure 1.Example of displaying the SA algorithm answer.

Figure 2 .
Figure 2. Rewrite the SA algorithm answer.

Figure 6 .
Figure 6.Signal-to-noise diagram of SA algorithm.

Figure 7 .
Figure 7. Time graph of problem-solving by GAMS and SA algorithm.

Figure 10 .
Figure 10.Sensitivity analysis of the objective function to the customers demand parameter.

Figure 11 .
Figure 11.Sensitivity analysis of the objective function to the payload capacity parameter.

Figure 12 .
Figure 12.Sensitivity analysis of battery capacity parameter.

Figure
Figure analysis of the cost of energy consumption parameter.

Figure 14 .
Figure 14.Sensitivity analysis of all parameters in this section.

Table 1 .
Notations for the mathematical model.Set of customers, charging stations, and depot node  + 1;  ′ 0 =  ′ ∪ { + 1}  ′ 0 ,  + 1 Set of customers, charging stations, and depot node  + 1;  ′ 0 =  ′ ∪ { + 1}  Types of electric taxis (economy, luxury, van) Parameters  Traveling distance from node  to node ; ∀,  ′ 0   The company's income from the unit of distance traveled by an electric taxi of the type   Unit cost of energy consumption   The maximum payload that the electric taxi of type  can carry   Battery capacity of the electric taxi of type   Total cargo weight demanded by customer node , if  ∈  ′ 0 →    = 0  The curb weight of the electric taxi The average speed of electric taxi  from node  to node   Travel time from node  to node ; ∀,  ′ 0, +1  Latest time to start service allowed at node ; ∀ ′ 0, +1  Earliest time to start service allowed at node ; ∀ ′

Table 3 .
Parameter levels of SA algorithm.

Table 4 .
Optimal values of SA algorithm parameters.

Table 5 .
Analysis of demand parameter and payload capacity changes.

Table 6 .
Analysis of payload capacity changes.

Table 7 .
Analysis of demand changes.

Table 8 .
The results of solving sample problems by GAMS and SA algorithm.

Table 9 .
Some problems solved with SA algorithm in MATLAB.

Table 10 .
The convergence diagram of the SA algorithm in Table9problems.

Table 11 .
Other problems solved with SA algorithm in MATLAB.