# A Hybrid Simulated Annealing Heuristic for Multistage Heterogeneous Fleet Scheduling with Fleet Sizing Decisions

Authors: Bing Li; Xinyu Yang; Hua Xuan
Journal: Journal of Advanced Transportation
 (2019)

---

## Abstract

This paper deals with multistage heterogeneous fleet scheduling with fleet sizing decisions (MHFS-FSD). This MHFS-FSD attempts to integrate vehicles allocation and fleet sizing decisions considering the vehicle routing of multiple vehicle types. The problem is formulated as mixed integer programming model. The matrix formulation denoting vehicle allocation scheme is explored according to the characteristic of this problem. Generating vehicle allocation scheme with greedy heuristic procedure (VA-GHP) as initial solution of problem is presented. The USP-IVA method to update the initial solution generated by VA-GHP approach is developed. And then, incorporating VA-GHP and USP-IVA into simulated annealing algorithm, a novel heuristic called HSAH-GHP&IVA is proposed. Finally, some experiments are designed to test the proposed heuristic and the results show that the heuristic can generate reasonably good solutions in short CPU times.


---

## Body

## 1. Introduction

Freight transport companies frequently face the problem of managing a fleet of vehicles which circulate on networks, being dispatched loaded or repositioned empty between various freight terminals of the network. There is recognition of the importance of managing fleets in trucking industry, railroad industry, and shipping industry.Prior research has most commonly focused on a homogeneous fleet of vehicles. But due to the large amount of loads and complex transport routing conditions, many conflicts will arise between homogeneous fleet of vehicles and heterogeneity transport demands when freight transport schemes are carried out. Instead of homogeneous fleet with substitution of single vehicle types for different types of transport demands, setting up heterogeneous fleets of vehicles is an effective strategy to enhance the efficiency of vehicle utilisation.In addition, the fleet size decisions are significant to improve the freight carry capacity of fleet. The insufficient vehicles and oversize vehicles can be ineffective in supporting fleet management.The study on heterogeneous fleet of various vehicle types with fleet sizing decisions is regarded as a promising perspective for the improvement of freight transport scheduling systems. The combinatorial nature of the problem causes the model formulated to be very hard to solve. Therefore, the research focuses on the development of powerful methods that are able to obtain solutions in reasonable time.In this paper, we focus on scheduling a heterogeneous fleet of various vehicles over time to serve a set of loads with fleet sizing decisions considering the vehicle routing of multiple vehicle types. We name this type of problem as multistage heterogeneous fleet scheduling with fleet sizing decisions (MHFS-FSD). The objective is to maximise total profits over the entire planning horizon. We develop a mixed integer programming model and present a novel hybrid simulated annealing heuristic called HSAH-GHP&IVA.The remainder of this paper is organised as follows. Section2 of this paper discusses related earlier research efforts. Section 3 is devoted to the mathematical description of the MHFS-FSD. In Section 4, we describe the mathematical formulation. Section 5 explains the approach of generating vehicle allocation scheme with greedy heuristic procedure (VA-GHP) as initial solution of problem. We develop the updating solution procedure by improving vehicle allocation (USP-IVA) in Section 6. And then the hybrid simulated annealing heuristic called HSAH-GHP&IVA is proposed in Section 7. The computational experiments are described in Section 8, and the effectiveness of the proposed method is shown from the computational results. The last section concludes with a summary of current work and extensions.


## 2. Literature Review

### 2.1. State of the Article

In this section, we review the relevant literature about fleet scheduling problems. Literature review indicates that fleet scheduling problems can be divided into four groups, that is, vehicle routing problem, vehicle assignment, transportation network optimisation, and fleet management.In recent years, many researches on vehicle routing problem have been carried out. The inventory routing problems with pickups and deliveries (IRP-PD) is studied by Archetti et al. [1] with a branch-and-cut method. Ciancio et al. [2] transformed the Mixed Capacitated General Routing Problem with Time Windows (MCGRPTW) into an equivalent node routing problem over a directed graph and solved the equivalent problem by using a branch-price-and-cut algorithm. A new combinatorial algorithm named OVRP_GELS based on gravitational emulation local search algorithm was obtained to solve the open vehicle routing problem (OVRP) by Hosseinabadi et al. [3] The Pickup and Delivery Traveling Salesman Problem with Multiple Stacks was researched by Sampaio and Urrutia [4]. They provided a new integer programming formulation with a polyhedral representation and a branch-and-cut algorithm for solving the proposed formulation. A two-phase approach to deal with the vehicle routing problems with backhauls and time windows is introduced by Reil et al. [5] They proposed Tabu search and first a multistart evolutionary strategy to minimise the total travel distance. Archetti et al. [6] designed and implemented a multistart heuristic which produces solutions with small errors when compared with optimal solutions obtained by solving an integer programming formulation with a commercial solver. Andelmin and Bartolini [7] presented an exact algorithm for solving the green vehicle routing problem and modeled the G-VRP as a set partitioning problem.The proper configuration of the vehicle assignment is the crucial point of the research on the optimisation of fleet resource allocation. Choi et al. [8] proposed a Dantzig-Wolfe decomposition approach for the vehicle assignment problem with demand uncertainty in a hybrid hub-and-spoke network. Spliet et al. [9] developed a branch-price-and-cut algorithm to arrange the time windows for time window assignment vehicle routing problems with the objective of minimising the expected transportation costs. Lin et al. [10] addressed dynamic vehicle allocation problem. Taking the waiting cost of transportation operation as the objective function of the problem, a Markov decision model is constructed and applied to the automated material handling system in semiconductor manufacturing. João et al. [11] presented a flight scheduling and fleet assignment optimisation model that may assist public authorities to establish the level of service requirements for subsidized air transport networks. Sargut et al. [12] introduced a multiobjective integrated crew rostering and vehicle assignment problem and developed a new multiobjective Tabu search algorithm. Li and Xuan [13] provided a greedy algorithm for solving the vehicle assignment with time window restraints.The optimisation of transportation network plays a vital role in improving fleet income and reducing various costs. The urban rail transit systems are optimised by Ozturk and Patrick [14]. They extended the setting to several stations and developed a heuristic method, two mixed integer models, and a constraint after proposing an approximation algorithm and a pseudopolynomial dynamic programming algorithm between two sites. Zhang and Zhang [15] explored the dynamic shortest path from a single source to a destination in a given traffic network. To optimise the journey time for the traveler, a novel dynamic shortest path algorithm based on hybridizing genetic and ant colony algorithms was developed. A feasible flow-based iterative algorithm named THTMTP-A is presented to deal with the two-level hierarchical time minimisation transportation problem by Xie et al. [16]. Suhng and Lee [17] developed a new Link-Based Single Tree Building Algorithm in order to reduce the slow execution speed problem of the multitree building algorithm for shortest path searching in an Urban Road Transportation Network and proved its usefulness by comparing the proposed one with other algorithms. The suitability of three different global optimisation methods which included the branch and reduce method, the branch and cut method, and the combination of global and local search strategies for specifically the exact optimum solution of the nonlinear transportation problem is studied by Klansek and Psunder [18].The fleet management is one of the most important as it is a major fixed investment for starting any business. Hashemi and Sattarvand [19] studied the different management systems of the open pit mining equipment including nondispatching, dispatching, and blending solutions for the Sungun copper mine. A dispatching simulation model with the objective function of minimising truck waiting times had been developed. Zhu et al. [20] focused on solving the scheduled service network design problem for freight rail transportation and proposed a heuristic solution methodology integrating slope scaling, a dynamic block-generation mechanism, long-term-memory-based perturbation strategies, and ellipsoidal search. Li et al. [21, 22], respectively, developed an alternating solution strategy for the stochastic dynamic fleet scheduling problem with variable period and presented a piecewise method by updating preset control parameters for dynamic working vehicle scheduling problem. Tierney et al. [23] studied a central problem in the liner shipping industry called the liner shipping fleet repositioning problem and introduced a simulated annealing algorithm for above problem. Hugo et al. [24] designed an approximate dynamic programming algorithm for large-scale fleet management.


### 2.2. Contributions

Overall, the objectives of this research are twofold. The first objective is to develop a mathematical model of MHFS-FSD. In the proposed model, the matrix formulation denoting vehicle allocation scheme is explored according to the characteristic of this problem. The second objective is to propose an efficient methodology for solving the model. Specifically, the contributions of this paper are as follows.(1)
We develop the mixed integer programming model (MIP). This model discusses the integration of vehicles allocation and fleet sizing decisions considering the vehicle routing of multiple vehicle types. Particular vehicle type is assigned to some particular transport routing according to the classification of vehicles. The matching of vehicle type and transport route is given. According to the characteristic of this problem, the matrix formulation denoting vehicle allocation scheme is provided.(2)
According to the specific structure of the mixed integer programming model, the approach of generating vehicle allocation scheme with greedy heuristic procedure (VA-GHP) as initial solution of problem is presented. On the basis of the initial solution generated by VA-GHP approach, the USP-IVA method to update the initial solution is developed, and then incorporating VA-GHP and USP-IVA into simulated annealing algorithm, a novel heuristic called HSAH-GHP&IVA is proposed.(3)
To evaluate the performance of the heuristics proposed, we generated randomly three sets of problem instances with different size of freight terminal and time periods, considering different number of vehicles types. Three sets of experiments are oriented to evaluate the performances of the hybrid simulated annealing heuristic. The results show that the proposed heuristic is able to obtain reasonably good solutions in short CPU times.


## 2.1. State of the Article

In this section, we review the relevant literature about fleet scheduling problems. Literature review indicates that fleet scheduling problems can be divided into four groups, that is, vehicle routing problem, vehicle assignment, transportation network optimisation, and fleet management.In recent years, many researches on vehicle routing problem have been carried out. The inventory routing problems with pickups and deliveries (IRP-PD) is studied by Archetti et al. [1] with a branch-and-cut method. Ciancio et al. [2] transformed the Mixed Capacitated General Routing Problem with Time Windows (MCGRPTW) into an equivalent node routing problem over a directed graph and solved the equivalent problem by using a branch-price-and-cut algorithm. A new combinatorial algorithm named OVRP_GELS based on gravitational emulation local search algorithm was obtained to solve the open vehicle routing problem (OVRP) by Hosseinabadi et al. [3] The Pickup and Delivery Traveling Salesman Problem with Multiple Stacks was researched by Sampaio and Urrutia [4]. They provided a new integer programming formulation with a polyhedral representation and a branch-and-cut algorithm for solving the proposed formulation. A two-phase approach to deal with the vehicle routing problems with backhauls and time windows is introduced by Reil et al. [5] They proposed Tabu search and first a multistart evolutionary strategy to minimise the total travel distance. Archetti et al. [6] designed and implemented a multistart heuristic which produces solutions with small errors when compared with optimal solutions obtained by solving an integer programming formulation with a commercial solver. Andelmin and Bartolini [7] presented an exact algorithm for solving the green vehicle routing problem and modeled the G-VRP as a set partitioning problem.The proper configuration of the vehicle assignment is the crucial point of the research on the optimisation of fleet resource allocation. Choi et al. [8] proposed a Dantzig-Wolfe decomposition approach for the vehicle assignment problem with demand uncertainty in a hybrid hub-and-spoke network. Spliet et al. [9] developed a branch-price-and-cut algorithm to arrange the time windows for time window assignment vehicle routing problems with the objective of minimising the expected transportation costs. Lin et al. [10] addressed dynamic vehicle allocation problem. Taking the waiting cost of transportation operation as the objective function of the problem, a Markov decision model is constructed and applied to the automated material handling system in semiconductor manufacturing. João et al. [11] presented a flight scheduling and fleet assignment optimisation model that may assist public authorities to establish the level of service requirements for subsidized air transport networks. Sargut et al. [12] introduced a multiobjective integrated crew rostering and vehicle assignment problem and developed a new multiobjective Tabu search algorithm. Li and Xuan [13] provided a greedy algorithm for solving the vehicle assignment with time window restraints.The optimisation of transportation network plays a vital role in improving fleet income and reducing various costs. The urban rail transit systems are optimised by Ozturk and Patrick [14]. They extended the setting to several stations and developed a heuristic method, two mixed integer models, and a constraint after proposing an approximation algorithm and a pseudopolynomial dynamic programming algorithm between two sites. Zhang and Zhang [15] explored the dynamic shortest path from a single source to a destination in a given traffic network. To optimise the journey time for the traveler, a novel dynamic shortest path algorithm based on hybridizing genetic and ant colony algorithms was developed. A feasible flow-based iterative algorithm named THTMTP-A is presented to deal with the two-level hierarchical time minimisation transportation problem by Xie et al. [16]. Suhng and Lee [17] developed a new Link-Based Single Tree Building Algorithm in order to reduce the slow execution speed problem of the multitree building algorithm for shortest path searching in an Urban Road Transportation Network and proved its usefulness by comparing the proposed one with other algorithms. The suitability of three different global optimisation methods which included the branch and reduce method, the branch and cut method, and the combination of global and local search strategies for specifically the exact optimum solution of the nonlinear transportation problem is studied by Klansek and Psunder [18].The fleet management is one of the most important as it is a major fixed investment for starting any business. Hashemi and Sattarvand [19] studied the different management systems of the open pit mining equipment including nondispatching, dispatching, and blending solutions for the Sungun copper mine. A dispatching simulation model with the objective function of minimising truck waiting times had been developed. Zhu et al. [20] focused on solving the scheduled service network design problem for freight rail transportation and proposed a heuristic solution methodology integrating slope scaling, a dynamic block-generation mechanism, long-term-memory-based perturbation strategies, and ellipsoidal search. Li et al. [21, 22], respectively, developed an alternating solution strategy for the stochastic dynamic fleet scheduling problem with variable period and presented a piecewise method by updating preset control parameters for dynamic working vehicle scheduling problem. Tierney et al. [23] studied a central problem in the liner shipping industry called the liner shipping fleet repositioning problem and introduced a simulated annealing algorithm for above problem. Hugo et al. [24] designed an approximate dynamic programming algorithm for large-scale fleet management.


## 2.2. Contributions

Overall, the objectives of this research are twofold. The first objective is to develop a mathematical model of MHFS-FSD. In the proposed model, the matrix formulation denoting vehicle allocation scheme is explored according to the characteristic of this problem. The second objective is to propose an efficient methodology for solving the model. Specifically, the contributions of this paper are as follows.(1)
We develop the mixed integer programming model (MIP). This model discusses the integration of vehicles allocation and fleet sizing decisions considering the vehicle routing of multiple vehicle types. Particular vehicle type is assigned to some particular transport routing according to the classification of vehicles. The matching of vehicle type and transport route is given. According to the characteristic of this problem, the matrix formulation denoting vehicle allocation scheme is provided.(2)
According to the specific structure of the mixed integer programming model, the approach of generating vehicle allocation scheme with greedy heuristic procedure (VA-GHP) as initial solution of problem is presented. On the basis of the initial solution generated by VA-GHP approach, the USP-IVA method to update the initial solution is developed, and then incorporating VA-GHP and USP-IVA into simulated annealing algorithm, a novel heuristic called HSAH-GHP&IVA is proposed.(3)
To evaluate the performance of the heuristics proposed, we generated randomly three sets of problem instances with different size of freight terminal and time periods, considering different number of vehicles types. Three sets of experiments are oriented to evaluate the performances of the hybrid simulated annealing heuristic. The results show that the proposed heuristic is able to obtain reasonably good solutions in short CPU times.


## 3. Problem Description and Analysis

This section describes the problem of scheduling a heterogeneous fleet of various vehicles over time to serve a set of loads with fleet sizing decisions. We make a brief overview of some foundational concepts and analysis of the problem.


### 3.1. Problem Description

The freight transportation scheduling frequently aims to obtain appropriate capacity of transportation systems, discover the potential of fleets of vehicles which circulate on networks, and improve transportation efficiency by arranging daily transportation production reasonably, organising vehicles flow adjustment scientifically, and making plans for loading and unloading operations in their destination terminal of the networks efficiently.This paper discusses the integration of vehicles allocation and fleet sizing decisions considering the vehicle routing of multiple vehicle types. LetG(N,E) represent the freight transport network, where N is a set of freight terminal andE represents the set of the movement of vehicles between a pair of freight terminal. We assume that time is divided into a set of discrete time periods H=h∣h=1,...,K whereK is the length of the planning horizon and h is the time period. All transitions require multiperiod travel times. Let tij be the travel time from i to j,i,j∈N. The travel time between any pair of terminals is different length of time period and integer multiple of one time period. The operators of transportation systems manage a heterogeneous fleet of various vehicles. We define the set of classification of vehicles as W=1,⋯,P. Particular vehicle type is assigned to some particular vehicle routing according to the classification of vehicles. The matching of vehicle type and vehicle routing is given. The use of operation vehicles will generate fixed expenses such as vehicle leasing and daily maintenance. It is crucial to decide the size of various vehicle types reasonably.To maximise total profits over a given horizon, the problem is to optimise the fleet sizing and make the vehicle utilisation decision in view of the profit derived from moving a load from terminali to terminaljdenoted as αijw, the cost of moving an empty vehicle from terminal i to terminal j denoted as βijw, the fixed costs of owning or leasing a type of vehicles denoted as γw, and vehicle routing of vehicle types denoted as πw. Figure 1 displays a snapshot of dynamic space-time network for fleet scheduling with multiperiod travel times.Figure 1
Dynamic space-time network for fleet scheduling with multiperiod travel times.


### 3.2. Understanding Vehicle Allocation between Pairs of Terminals

A fleet of vehicles must be assigned to loads that have the effect of moving the vehicle from one terminal to the next. The operator of the fleet must either assign each vehicle to a requested loaded movement or move it empty to another terminal to pick up a requested loaded movement, or to hold it in the same terminal until the next time period. So the allocation problem of loaded vehicle and empty vehicle arises.Serving loads results in the relocation of vehicles. Demands request a vehicle to be available in a specific terminal on a specific time period to carry the load to a given destination terminal. Each load must be served by a vehicle. Loaded movement will generate revenue. So the loaded vehicle assignment is the basis of the vehicle allocation.The load for movements between various terminals is often imbalanced, and this implies the need for redistribution of empty vehicles over the transport network from terminals at which they have become idle to terminals at which they can be reused. It is important to manage empty vehicle flows in transport network. The loads create requirements for empty vehicle. The arrivals of loaded vehicle flow create the supplies of empties available for repositioning.Vehicle allocation has significant effects on the vehicle distribution or utilisation decision and then effects on capacity and efficiency of transport network.


### 3.3. Understanding Heterogeneous Fleet of Vehicles

Due to the huge transport networks, the large amount of loads, and complex transport route conditions, many conflicts will arise between homogeneous fleet of vehicles and heterogeneity transport demands when freight transport schemes are carried out. A homogeneous fleet with single vehicle type means inefficiency in satisfying the diverse transport demand. When the operators manage a homogeneous fleet, the vehicle shortage and idleness may occur synchronously. The worsening trend of vehicle utilisation is also possible to appear so that transport efficiency is reduced and transport cost is enhanced. Instead of homogeneous fleet with substitution of single vehicle types for different types of transport demands, setting up heterogeneous fleets of vehicles is an effective strategy to eliminate the above unreasonable phenomena, optimise the vehicle allocation, and enhance the efficiency of vehicle utilisation.Meanwhile, the fleet size decisions are significant to improve the freight carry capacity of fleet. The insufficient vehicles seriously affect the implementation of transport work. If the number of available vehicles is not enough to meet transport demand, the penalty cost for unmet demand will generate. On the other hand, the oversized supply of fleet far more than the transport demand causes a large number of vehicles to be idle.


### 3.4. Understanding Vehicle Routing for Particular Vehicle Types

Many works have assumed a homogeneous fleet of vehicles. But most problems in practice have multiple vehicle types with limited substitution. In this paper, we extend the methodology to handle heterogeneous fleets of vehicles with different vehicle types for different transport routing.The two stage freight transport network is comprised of freight supply location as origin terminal and freight demand location as destination terminal. In practical considerations, the route from any supply location to demand location cannot always be used to carry freight due to limited vehicle types, transport routing condition, transport routing external environment, and so on. The specific reasons can be diverse as follows:(1) Some transport routes cannot be used by some specific vehicle type due to carry capacity of the road.(2) Some specific vehicle types are not allowed to travel on certain transport routes due to physical condition of the road versus maintenance condition of vehicles.(3) Due to inexistence of supply location and demand location for freight exchange, some transport routes disappear.An illustration of the vehicle routing of multiple vehicle types is shown in Figure2.Figure 2
The vehicle routing of multiple vehicle types.Because some specific vehicle types are only allowed to travel on certain transport routes, we useπ as the matching matrix such that πijw=1 if w types of vehicle can be assigned to movement from terminal i to j, and πijw=0 otherwise. We present the notation for more clearly expressing the matching matrix here.N is the set of freight terminal in the freight transport network. Define N=1,⋯,A whereA is the number of freight terminal.W is the set of vehicle types in the heterogeneous fleet. Define W=1,⋯,P where P is the total number of available vehicles of various types.T is the matrix of travel time of vehicles set between a pair of freight terminal. We have(1)T=0⋯t1j⋯t1A⋮⋮⋮ti1⋯tij⋯tiA⋮⋮⋮tA1⋯tAj⋯0,∀i,j∈Nπ is the matrix of vehicle routing for certain vehicle type. The matrix can be expressed by π=πw∣w=1,⋯,P. The parameter πw is used to represent practical situations where vehicle type w can be allowed to transport load on certain routes. Correspondingly, we denote that(2)πw=π11w⋯π1jw⋯π1Aw⋮⋮⋮πi1w⋯πijw⋯πiAw⋮⋮⋮πA1w⋯πAjw⋯πAAwπijw=0,if  vehicle  typewcan  move  from  terminalitoj1,otherwise∀i,j∈N,w∈W.


## 3.1. Problem Description

The freight transportation scheduling frequently aims to obtain appropriate capacity of transportation systems, discover the potential of fleets of vehicles which circulate on networks, and improve transportation efficiency by arranging daily transportation production reasonably, organising vehicles flow adjustment scientifically, and making plans for loading and unloading operations in their destination terminal of the networks efficiently.This paper discusses the integration of vehicles allocation and fleet sizing decisions considering the vehicle routing of multiple vehicle types. LetG(N,E) represent the freight transport network, where N is a set of freight terminal andE represents the set of the movement of vehicles between a pair of freight terminal. We assume that time is divided into a set of discrete time periods H=h∣h=1,...,K whereK is the length of the planning horizon and h is the time period. All transitions require multiperiod travel times. Let tij be the travel time from i to j,i,j∈N. The travel time between any pair of terminals is different length of time period and integer multiple of one time period. The operators of transportation systems manage a heterogeneous fleet of various vehicles. We define the set of classification of vehicles as W=1,⋯,P. Particular vehicle type is assigned to some particular vehicle routing according to the classification of vehicles. The matching of vehicle type and vehicle routing is given. The use of operation vehicles will generate fixed expenses such as vehicle leasing and daily maintenance. It is crucial to decide the size of various vehicle types reasonably.To maximise total profits over a given horizon, the problem is to optimise the fleet sizing and make the vehicle utilisation decision in view of the profit derived from moving a load from terminali to terminaljdenoted as αijw, the cost of moving an empty vehicle from terminal i to terminal j denoted as βijw, the fixed costs of owning or leasing a type of vehicles denoted as γw, and vehicle routing of vehicle types denoted as πw. Figure 1 displays a snapshot of dynamic space-time network for fleet scheduling with multiperiod travel times.Figure 1
Dynamic space-time network for fleet scheduling with multiperiod travel times.


## 3.2. Understanding Vehicle Allocation between Pairs of Terminals

A fleet of vehicles must be assigned to loads that have the effect of moving the vehicle from one terminal to the next. The operator of the fleet must either assign each vehicle to a requested loaded movement or move it empty to another terminal to pick up a requested loaded movement, or to hold it in the same terminal until the next time period. So the allocation problem of loaded vehicle and empty vehicle arises.Serving loads results in the relocation of vehicles. Demands request a vehicle to be available in a specific terminal on a specific time period to carry the load to a given destination terminal. Each load must be served by a vehicle. Loaded movement will generate revenue. So the loaded vehicle assignment is the basis of the vehicle allocation.The load for movements between various terminals is often imbalanced, and this implies the need for redistribution of empty vehicles over the transport network from terminals at which they have become idle to terminals at which they can be reused. It is important to manage empty vehicle flows in transport network. The loads create requirements for empty vehicle. The arrivals of loaded vehicle flow create the supplies of empties available for repositioning.Vehicle allocation has significant effects on the vehicle distribution or utilisation decision and then effects on capacity and efficiency of transport network.


## 3.3. Understanding Heterogeneous Fleet of Vehicles

Due to the huge transport networks, the large amount of loads, and complex transport route conditions, many conflicts will arise between homogeneous fleet of vehicles and heterogeneity transport demands when freight transport schemes are carried out. A homogeneous fleet with single vehicle type means inefficiency in satisfying the diverse transport demand. When the operators manage a homogeneous fleet, the vehicle shortage and idleness may occur synchronously. The worsening trend of vehicle utilisation is also possible to appear so that transport efficiency is reduced and transport cost is enhanced. Instead of homogeneous fleet with substitution of single vehicle types for different types of transport demands, setting up heterogeneous fleets of vehicles is an effective strategy to eliminate the above unreasonable phenomena, optimise the vehicle allocation, and enhance the efficiency of vehicle utilisation.Meanwhile, the fleet size decisions are significant to improve the freight carry capacity of fleet. The insufficient vehicles seriously affect the implementation of transport work. If the number of available vehicles is not enough to meet transport demand, the penalty cost for unmet demand will generate. On the other hand, the oversized supply of fleet far more than the transport demand causes a large number of vehicles to be idle.


## 3.4. Understanding Vehicle Routing for Particular Vehicle Types

Many works have assumed a homogeneous fleet of vehicles. But most problems in practice have multiple vehicle types with limited substitution. In this paper, we extend the methodology to handle heterogeneous fleets of vehicles with different vehicle types for different transport routing.The two stage freight transport network is comprised of freight supply location as origin terminal and freight demand location as destination terminal. In practical considerations, the route from any supply location to demand location cannot always be used to carry freight due to limited vehicle types, transport routing condition, transport routing external environment, and so on. The specific reasons can be diverse as follows:(1) Some transport routes cannot be used by some specific vehicle type due to carry capacity of the road.(2) Some specific vehicle types are not allowed to travel on certain transport routes due to physical condition of the road versus maintenance condition of vehicles.(3) Due to inexistence of supply location and demand location for freight exchange, some transport routes disappear.An illustration of the vehicle routing of multiple vehicle types is shown in Figure2.Figure 2
The vehicle routing of multiple vehicle types.Because some specific vehicle types are only allowed to travel on certain transport routes, we useπ as the matching matrix such that πijw=1 if w types of vehicle can be assigned to movement from terminal i to j, and πijw=0 otherwise. We present the notation for more clearly expressing the matching matrix here.N is the set of freight terminal in the freight transport network. Define N=1,⋯,A whereA is the number of freight terminal.W is the set of vehicle types in the heterogeneous fleet. Define W=1,⋯,P where P is the total number of available vehicles of various types.T is the matrix of travel time of vehicles set between a pair of freight terminal. We have(1)T=0⋯t1j⋯t1A⋮⋮⋮ti1⋯tij⋯tiA⋮⋮⋮tA1⋯tAj⋯0,∀i,j∈Nπ is the matrix of vehicle routing for certain vehicle type. The matrix can be expressed by π=πw∣w=1,⋯,P. The parameter πw is used to represent practical situations where vehicle type w can be allowed to transport load on certain routes. Correspondingly, we denote that(2)πw=π11w⋯π1jw⋯π1Aw⋮⋮⋮πi1w⋯πijw⋯πiAw⋮⋮⋮πA1w⋯πAjw⋯πAAwπijw=0,if  vehicle  typewcan  move  from  terminalitoj1,otherwise∀i,j∈N,w∈W.


## 4. Mathematical Formulation

In this section, we define the notation and formulate multistage heterogeneous fleet scheduling with fleet sizing decisions (MHFS-FSD) as a mixed integer programming model. In the following we describe the model parameters and variables used to model the MHFS-FSD.


### 4.1. Model Parameters

The parameters are needed to describe the problem.t i j is the travel time moving from origin i to destination j and is the integer multiple of the scheduling time period.π i j w is the parameter for representing whether or not vehicle type w can be allowed to carry freight along transport routing (i,j). The value of πijw is 1 if the w types of vehicle can carry freight from the terminal i to the terminal j; otherwise, the value is 0, ∀w∈W, and ∀i,j∈N.q i j h is the set of loads with origin i and destination j having h as a feasible departure, ∀h∈H, ∀i,j∈N.α i j is the profits generated by a unit load being dispatched from terminal i to terminal j, ∀w∈W, ∀i,j∈N.β i j w is the cost of w type vehicle being repositioned empty from terminal i to terminal j, ∀w∈W, ∀i,j∈N.γ w is the purchase or lease cost for the w type of vehicle, ∀w∈W.λ w is the loading capacity of w type vehicle ∀w∈W. That is to say, it is the amount of loads which can be serviced by one w type vehicle.z i t w is the number of available vehicle of w type at terminal i in time periodt, ∀w∈W, ∀i∈N.


### 4.2. Model Variables

The decision variables of the mathematical model arex i j h , w is the number of w type of loaded vehicles from terminal i to terminal j at time period h, ∀i,j∈N, h∈H.y i j h , w is the number of w type of empty vehicles from terminal i to terminal j at time period h, ∀i,j∈N,h∈H.z w is the number for available vehicle of w type.


### 4.3. Mathematical Models

The mathematical models for multistage heterogeneous fleet scheduling with fleet sizing decisions (MHFS-FSD) are given by(3)maxFx,y,z=∑w∈W∑h=1K∑i∈N∑j∈Nαij·min⁡qijh,∑w∈Wλwxijh,w-βijwxijh,w-βijwyijh,w-γwzw(4)s.t.∑j∈Nxijh,w+yijh,wπijw=zih,w,∀i∈N,∀h∈H,∀w∈W(5)∑w∈Wλwxijh,w≤qijh,∀i,j∈N,∀h∈H(6)∑i∈Nxijh-tij,w+yijh-tij,wπijw=zjh,w,∀j∈N,∀h∈H,∀w∈W(7)∑i∈Nzih,w=zw,∀h∈H(8)xijh,w≥0,yijh,w≥0,zw≥0,integer,∀i,j∈N,∀h∈H,∀w∈WThe MHFS-FSD model belongs to a mixed integer programming model (MIP). The objective function (3) aims to maximise the total profit generated by the fleet scheduling scheme during the whole planning horizon. Constraint (4) denotes the number of dispatching loaded vehicles and empty vehicles are determined by supplies of vehicles in time period. Constraint (5) ensures that the number of loaded vehicle cannot exceed available requirements of all loads in each terminal at each time period. Constraint (6) denotes updating equation for total supply of vehicle of various types. Constraint (7) shows having flow conservation of vehicle in each type. Finally, constraint (8) defines the range of decision variables constraints.As an NP problem, using traditional method for solving MHFS-FSD model is difficult and inefficient. A novel heuristic method is considered to solve this problem.


## 4.1. Model Parameters

The parameters are needed to describe the problem.t i j is the travel time moving from origin i to destination j and is the integer multiple of the scheduling time period.π i j w is the parameter for representing whether or not vehicle type w can be allowed to carry freight along transport routing (i,j). The value of πijw is 1 if the w types of vehicle can carry freight from the terminal i to the terminal j; otherwise, the value is 0, ∀w∈W, and ∀i,j∈N.q i j h is the set of loads with origin i and destination j having h as a feasible departure, ∀h∈H, ∀i,j∈N.α i j is the profits generated by a unit load being dispatched from terminal i to terminal j, ∀w∈W, ∀i,j∈N.β i j w is the cost of w type vehicle being repositioned empty from terminal i to terminal j, ∀w∈W, ∀i,j∈N.γ w is the purchase or lease cost for the w type of vehicle, ∀w∈W.λ w is the loading capacity of w type vehicle ∀w∈W. That is to say, it is the amount of loads which can be serviced by one w type vehicle.z i t w is the number of available vehicle of w type at terminal i in time periodt, ∀w∈W, ∀i∈N.


## 4.2. Model Variables

The decision variables of the mathematical model arex i j h , w is the number of w type of loaded vehicles from terminal i to terminal j at time period h, ∀i,j∈N, h∈H.y i j h , w is the number of w type of empty vehicles from terminal i to terminal j at time period h, ∀i,j∈N,h∈H.z w is the number for available vehicle of w type.


## 4.3. Mathematical Models

The mathematical models for multistage heterogeneous fleet scheduling with fleet sizing decisions (MHFS-FSD) are given by(3)maxFx,y,z=∑w∈W∑h=1K∑i∈N∑j∈Nαij·min⁡qijh,∑w∈Wλwxijh,w-βijwxijh,w-βijwyijh,w-γwzw(4)s.t.∑j∈Nxijh,w+yijh,wπijw=zih,w,∀i∈N,∀h∈H,∀w∈W(5)∑w∈Wλwxijh,w≤qijh,∀i,j∈N,∀h∈H(6)∑i∈Nxijh-tij,w+yijh-tij,wπijw=zjh,w,∀j∈N,∀h∈H,∀w∈W(7)∑i∈Nzih,w=zw,∀h∈H(8)xijh,w≥0,yijh,w≥0,zw≥0,integer,∀i,j∈N,∀h∈H,∀w∈WThe MHFS-FSD model belongs to a mixed integer programming model (MIP). The objective function (3) aims to maximise the total profit generated by the fleet scheduling scheme during the whole planning horizon. Constraint (4) denotes the number of dispatching loaded vehicles and empty vehicles are determined by supplies of vehicles in time period. Constraint (5) ensures that the number of loaded vehicle cannot exceed available requirements of all loads in each terminal at each time period. Constraint (6) denotes updating equation for total supply of vehicle of various types. Constraint (7) shows having flow conservation of vehicle in each type. Finally, constraint (8) defines the range of decision variables constraints.As an NP problem, using traditional method for solving MHFS-FSD model is difficult and inefficient. A novel heuristic method is considered to solve this problem.


## 5. Generating Initial Solution with VA-GHP

According to the characteristic of this problem, the matrix formulation denoting vehicle allocation scheme is provided. And then the approach of generating vehicle allocation scheme with greedy heuristic procedure as initial solution of problem is presented.


### 5.1. Matrix Formulation of Vehicle Allocation

The solution of MHFS-FSD model represents the vehicle allocation scheme. The solution is expressed as two stage structure. The first term in the solution is matrix formulation of loaded vehicle dispatching scheme. The second term is the matrix formulation of empty vehicle reposition scheme.(  1) The Matrix Formulation of Loaded Vehicle Dispatching Scheme. The first section of the solution is used to represent the loaded vehicle dispatching scheme of carrying freight from the origin terminal to the destination terminal. This section is formulated as the matrix form which has A rows andA columns. Each component in row and row intersection of the matrix indicates the number of w type loaded vehicles from origin terminal to destination terminal. The allocation scheme of loaded vehicle is written as the following formula:(9)xw=x11w⋯x1jw⋯x1Aw⋮⋮⋮xi1w⋯xijw⋯xiAw⋮⋮⋮xA1w⋯xAjw⋯xAAw,w∈W.(  2) The Matrix Formulation of Empty Vehicle Reposition Scheme. The second section of the solution is used to denote the empty vehicle reposition scheme of vehicle moving empty from the origin terminal to the destination terminal. This section is also formulated as the matrix form which has A rows and A columns. Each component in row and row intersection of the matrix expresses the number of w type empty vehicles from origin terminal to destination terminal. The allocation scheme of loaded vehicle is expressed as the following formula:(10)yw=y11w⋯y1jw⋯y1Aw⋮⋮⋮yi1w⋯yijw⋯yiAw⋮⋮⋮yA1w⋯yAjw⋯yAAw,w∈W


### 5.2. Greedy Heuristic Procedure to Generate Vehicle Allocation Scheme

According to the fundamental principles of prioritising to consider the vehicle types with large load capacity and transport routing with large carry capacity, the greedy heuristic procedure is presented to generate vehicle allocation scheme. We name this approach as vehicle allocation with greedy heuristic procedure (VA-GHP). An overview of the framework of VA-GHP is explained as follows.Stage 1 (set of vehicle types selection sequence).
According to the vehicle distributionzih,w in terminal i at h time period, the set of vehicle types selection sequence Wih is obtained with ranking the vehicle types in descending order by loading capacity λw. We denote δ as serial number of the vehicle type in this set and δ^ as the total quantity of available vehicle types. Thus, Wih has the form of(11)Wih=wihδ∣wihδ=w,w∈W,δ=1,⋯,δ^Stage 2 (set of fleet size selection sequence).
LetZih be the set of fleet size selection sequence. According to the set of vehicle types selection sequence Wih, we can express Zih as(12)Zih=zihδ∣zihδ=zih,w,w∈Wih,δ=1,⋯,δ^Stage 3 (set of transport routing selection sequence).
We randomly choose the terminali from the set of freight terminal and compute the revenue αijqijh created by loaded vehicle dispatched from origin i to any available destination j. We rank the transport routing in order of the revenue from high to low and the set of transport routing selection sequence Sih can be given. We use ξ to denote the serial number of transport routing in this set and the total number of transport routing is represented as ξ^ih. Here Sih is written as(13)Sih=sihξ∣sihξ=i,j,j∈N,ξ=1,⋯,ξ^ihThe set of transport routing in terminali ath time period is described in Figure 3.Figure 3
A diagram of the set of transport routing selection sequence.Stage 4 (loaded vehicles dispatching scheme).
⁡
S tep 4.1(loaded vehicles dispatching scheme of transport routing). The transport routing sih(ξ) is selected in arranged order from the set of transport routing selection sequence Sih. According to the vehicle routing matrix πw for certain vehicle type, the available vehicle types wih(δ) allowed to travel on the transport routing sih(ξ) are chosen from the vehicle type selection sequence Wih. And the fleet size zih(δ) of corresponding vehicle types is determined from fleet size selection sequence Zih. We use the formula (14) to compute loaded vehicles dispatching scheme:(14)xih,δξ=min⁡qihξλδ,zihδ
S tep 4.2(adjusting available loads number and available fleet size for transport routing). (1) If there are qih(ξ)/λδ≥zih(δ), then xih,δ(ξ)=zih(δ). Besides, let qih(ξ)=qih(ξ)-λδzih(δ) and zih(δ)=0; turn to Step4.3.
(2) If there areqih(ξ)/λδ≤zih(δ), then xih,δ(ξ)=qih(ξ)/λδ. Besides, let qih(ξ)=0 and zih(δ)=zih(δ)-qih(ξ)/λδ; turn to Step4.4.
S tep 4.3(updating three selection sequence sets).
⁡
(  1) Updating the Set of Transport Routing Selection Sequence. According to available updated loads number qih(ξ), the revenue of transport routing with an origin terminali is recalculated. The set of transport routing selection sequence S-ih is updated with ranking the transport routing in descending order of the revenue αijqih(ξ). The new set S-ih is expressed as follows:(15)S-ih=s-ihξ-∣s-ihξ-=i,j-,j-∈N,ξ-=1,⋯,ξ^ih
(  2) Updating the Set of Vehicle Types Selection Sequence. Removing the vehicle types wih(δ) that have been used up from vehicle selection sequence set Wih, the new set of vehicle type selection sequence W-ih can be given by(16)W-ih=Wih-wihδ=w-ihδ-∣w-ihδ-=wihδ,δ-=δ+1,⋯,δ^
(  3) Updating the Set of Fleet Size Selection Sequence. Removing the fleet size zih(δ) of the vehicle types wih(δ) from fleet size selection sequence set Zih, the new set of fleet size selection sequence Z-ih can be given by(17)Z-ih=Zih-zihδ=z-ihδ-∣z-ihδ-=zihδ,δ=δ+1,⋯,δ^
S tep 4.4(updating two selection sequence sets).
⁡
(  1) Updating the Set of Transport Routing Selection Sequence. The new set of transport routing selection sequence S-ih is given with removing the transport routing sih(ξ) from transport routing selection sequence Sih. The new set S-ih is expressed as follows:(18)S-ih=Sih-sihξ=s-ihξ-∣s-ihξ-=sihξ,ξ-=ξ+1,⋯,ξ^ih
(  2) Updating the Set of Fleet Sizing Selection Sequence. The new set of fleet size selection sequence Z-ih is obtained by updating the corresponding fleet size zih(δ) of vehicle type wih(δ) in the fleet size selection sequence set Zih. The new set Z-ih is denoted as(19)Z-ih=z-ihδ-∣z-ihδ=zihδ-qihξλδ,z-ihδ+1=zihδ+1,δ-=δ,⋯,δ^Stage 5 (empty vehicle reposition scheme).
We can obtain the loaded vehicles dispatching schemexij(h+1-tij),w by the procedure described in Stage 4. Substituting them into equation ∑i∈Nxij(h+1-tij),wπijw=z~jh+1,w, the vehicles distribution z~jh+1,w can be obtained at terminalj at time period h+1. And, according to the load number qijh+1 at time period h+1, empty vehicle reposition scheme yijh,w can be obtained.S tep 5.1 (calculating the number of loads waiting to be arranged). The number of loads waiting to be arranged at each freight terminal q~jh+1 is computed by the following formula:(20)q~jh+1=∑k∈Nqjkh+1-∑w∈Wλwz~jh+1,wS tep 5.2 (set of freight terminal selected sequence with lack of empty vehicle). For the freight terminal where the number of load waiting to be arranged q~jh+1 is nonpositive, it means that there are empty vehicles enough not to continue to call the empty vehicles. Following the abovementioned principles, the freight terminals with lack of empty vehicle are ranked in descending order of the revenues ∑k∈Nαjkq~jkh+1. The set of freight terminal selected sequence with lack of empty vehicle J~h+1 is formed. We denote σ as the serial number of freight terminal with lack of empty vehicle in this set J~h+1 and σ^ as the total quantity of freight terminal with lack of empty vehicle. So the new set J~h+1 can be expressed in the form by(21)J~h+1=jσ∣jσ=j,j∈N,σ=1⋯,σ^S tep 5.3 (set of empty vehicles shortage amount selected sequence). We count the shortage amount of freight terminal with lack of empty vehicle in the set J~h+1. Further, the set of empty vehicles shortage amount selected sequence Q~jh+1 is obtained and is denoted as Q~jh+1=q~jh+1,j∈J~.S tep 5.4 (set of freight terminal selected sequence for dispatching empty vehicle). The destination terminaljis successively taken out from the set of freight terminal selected sequence with lack of empty vehicle J~h+1. Then, according to the vehicle routing matrix πw, the origin terminal i that has a link with destination terminal j is taken out. And, ranking the origin terminal i in ascending order of the distance, the set of freight terminal selected sequence for dispatching empty vehicle I~jh is formed. We denote θ as the serial number of freight terminal for dispatching empty vehicle in this set I~jh and σ^ as the total quantity of freight terminal for dispatching empty vehicle. So the new set I~jh can be expressed in the form by(22)I~jh=iθ∣πiθjw=1,iθ∈N,w∈W,j∈J~h+1,θ=1,⋯,θ^jhS tep 5.5 (empty vehicle reposition scheme). We choose the destination terminaljfrom the set of freight terminal selected sequence with lack of empty vehicle J~h+1 and the origin terminal i(θ) from the corresponding set of freight terminal selected sequence for dispatching empty vehicle I~jh. In accordance with vehicle distribution Z-i(θ)h at freight terminal for dispatching empty vehicle i(θ) and loads number q~jh+1 at freight terminal with lack of empty vehicle, we make use of the principle of priority arranging large capacity vehicle types to conduct the empty vehicle reposition as the following formula:(23)yjh,δθ=min⁡q~jh+1λδ,z-iθhδ,i∈I~jh(1) If there areq~jh+1/λδ≥z-i(θ)h(δ), then yjh,δ(θ)=z-i(θ)h(δ). Additionally, let q~jh+1=q~jh+1-λδz-i(θ)h(δ) and z-i(θ)h(δ)=0. If there are still some other vehicle types unarranged at freight terminal for dispatching empty vehicle i(θ), they are furtherly made by using equation (23). Otherwise, succeeding freight terminal i(θ) should be taken out from the set of freight terminal selected sequence for dispatching empty vehicle I~jh.(2) If there areq~jh+1/λδ≤z-i(θ)h(δ), then yjh,δ(θ)=q~jh+1/λδ. Besides, let q~jh+1=0 and z-i(θ)h(δ)=z-i(θ)h(δ)-q~jh+1/λδ. In addition, next destination stationj is taken out from set of freight terminal selected sequence with lack of empty vehicle J~h+1.Stage 6 (updating the vehicle distributions).
We use∑i∈Nxijh+1-tij,w+yijh+1-tij,wπijw=zjh+1,w to count the vehicle distributions Zjh+1,w in h period and return to Stage1.The approach of vehicle allocation with greedy heuristic procedure (VA-GHP) is shown in Figure4.Figure 4
Loaded vehicles dispatching and empty vehicle reposition scheme with greedy heuristic procedure.


## 5.1. Matrix Formulation of Vehicle Allocation

The solution of MHFS-FSD model represents the vehicle allocation scheme. The solution is expressed as two stage structure. The first term in the solution is matrix formulation of loaded vehicle dispatching scheme. The second term is the matrix formulation of empty vehicle reposition scheme.(  1) The Matrix Formulation of Loaded Vehicle Dispatching Scheme. The first section of the solution is used to represent the loaded vehicle dispatching scheme of carrying freight from the origin terminal to the destination terminal. This section is formulated as the matrix form which has A rows andA columns. Each component in row and row intersection of the matrix indicates the number of w type loaded vehicles from origin terminal to destination terminal. The allocation scheme of loaded vehicle is written as the following formula:(9)xw=x11w⋯x1jw⋯x1Aw⋮⋮⋮xi1w⋯xijw⋯xiAw⋮⋮⋮xA1w⋯xAjw⋯xAAw,w∈W.(  2) The Matrix Formulation of Empty Vehicle Reposition Scheme. The second section of the solution is used to denote the empty vehicle reposition scheme of vehicle moving empty from the origin terminal to the destination terminal. This section is also formulated as the matrix form which has A rows and A columns. Each component in row and row intersection of the matrix expresses the number of w type empty vehicles from origin terminal to destination terminal. The allocation scheme of loaded vehicle is expressed as the following formula:(10)yw=y11w⋯y1jw⋯y1Aw⋮⋮⋮yi1w⋯yijw⋯yiAw⋮⋮⋮yA1w⋯yAjw⋯yAAw,w∈W


## 5.2. Greedy Heuristic Procedure to Generate Vehicle Allocation Scheme

According to the fundamental principles of prioritising to consider the vehicle types with large load capacity and transport routing with large carry capacity, the greedy heuristic procedure is presented to generate vehicle allocation scheme. We name this approach as vehicle allocation with greedy heuristic procedure (VA-GHP). An overview of the framework of VA-GHP is explained as follows.Stage 1 (set of vehicle types selection sequence).
According to the vehicle distributionzih,w in terminal i at h time period, the set of vehicle types selection sequence Wih is obtained with ranking the vehicle types in descending order by loading capacity λw. We denote δ as serial number of the vehicle type in this set and δ^ as the total quantity of available vehicle types. Thus, Wih has the form of(11)Wih=wihδ∣wihδ=w,w∈W,δ=1,⋯,δ^Stage 2 (set of fleet size selection sequence).
LetZih be the set of fleet size selection sequence. According to the set of vehicle types selection sequence Wih, we can express Zih as(12)Zih=zihδ∣zihδ=zih,w,w∈Wih,δ=1,⋯,δ^Stage 3 (set of transport routing selection sequence).
We randomly choose the terminali from the set of freight terminal and compute the revenue αijqijh created by loaded vehicle dispatched from origin i to any available destination j. We rank the transport routing in order of the revenue from high to low and the set of transport routing selection sequence Sih can be given. We use ξ to denote the serial number of transport routing in this set and the total number of transport routing is represented as ξ^ih. Here Sih is written as(13)Sih=sihξ∣sihξ=i,j,j∈N,ξ=1,⋯,ξ^ihThe set of transport routing in terminali ath time period is described in Figure 3.Figure 3
A diagram of the set of transport routing selection sequence.Stage 4 (loaded vehicles dispatching scheme).
⁡
S tep 4.1(loaded vehicles dispatching scheme of transport routing). The transport routing sih(ξ) is selected in arranged order from the set of transport routing selection sequence Sih. According to the vehicle routing matrix πw for certain vehicle type, the available vehicle types wih(δ) allowed to travel on the transport routing sih(ξ) are chosen from the vehicle type selection sequence Wih. And the fleet size zih(δ) of corresponding vehicle types is determined from fleet size selection sequence Zih. We use the formula (14) to compute loaded vehicles dispatching scheme:(14)xih,δξ=min⁡qihξλδ,zihδ
S tep 4.2(adjusting available loads number and available fleet size for transport routing). (1) If there are qih(ξ)/λδ≥zih(δ), then xih,δ(ξ)=zih(δ). Besides, let qih(ξ)=qih(ξ)-λδzih(δ) and zih(δ)=0; turn to Step4.3.
(2) If there areqih(ξ)/λδ≤zih(δ), then xih,δ(ξ)=qih(ξ)/λδ. Besides, let qih(ξ)=0 and zih(δ)=zih(δ)-qih(ξ)/λδ; turn to Step4.4.
S tep 4.3(updating three selection sequence sets).
⁡
(  1) Updating the Set of Transport Routing Selection Sequence. According to available updated loads number qih(ξ), the revenue of transport routing with an origin terminali is recalculated. The set of transport routing selection sequence S-ih is updated with ranking the transport routing in descending order of the revenue αijqih(ξ). The new set S-ih is expressed as follows:(15)S-ih=s-ihξ-∣s-ihξ-=i,j-,j-∈N,ξ-=1,⋯,ξ^ih
(  2) Updating the Set of Vehicle Types Selection Sequence. Removing the vehicle types wih(δ) that have been used up from vehicle selection sequence set Wih, the new set of vehicle type selection sequence W-ih can be given by(16)W-ih=Wih-wihδ=w-ihδ-∣w-ihδ-=wihδ,δ-=δ+1,⋯,δ^
(  3) Updating the Set of Fleet Size Selection Sequence. Removing the fleet size zih(δ) of the vehicle types wih(δ) from fleet size selection sequence set Zih, the new set of fleet size selection sequence Z-ih can be given by(17)Z-ih=Zih-zihδ=z-ihδ-∣z-ihδ-=zihδ,δ=δ+1,⋯,δ^
S tep 4.4(updating two selection sequence sets).
⁡
(  1) Updating the Set of Transport Routing Selection Sequence. The new set of transport routing selection sequence S-ih is given with removing the transport routing sih(ξ) from transport routing selection sequence Sih. The new set S-ih is expressed as follows:(18)S-ih=Sih-sihξ=s-ihξ-∣s-ihξ-=sihξ,ξ-=ξ+1,⋯,ξ^ih
(  2) Updating the Set of Fleet Sizing Selection Sequence. The new set of fleet size selection sequence Z-ih is obtained by updating the corresponding fleet size zih(δ) of vehicle type wih(δ) in the fleet size selection sequence set Zih. The new set Z-ih is denoted as(19)Z-ih=z-ihδ-∣z-ihδ=zihδ-qihξλδ,z-ihδ+1=zihδ+1,δ-=δ,⋯,δ^Stage 5 (empty vehicle reposition scheme).
We can obtain the loaded vehicles dispatching schemexij(h+1-tij),w by the procedure described in Stage 4. Substituting them into equation ∑i∈Nxij(h+1-tij),wπijw=z~jh+1,w, the vehicles distribution z~jh+1,w can be obtained at terminalj at time period h+1. And, according to the load number qijh+1 at time period h+1, empty vehicle reposition scheme yijh,w can be obtained.S tep 5.1 (calculating the number of loads waiting to be arranged). The number of loads waiting to be arranged at each freight terminal q~jh+1 is computed by the following formula:(20)q~jh+1=∑k∈Nqjkh+1-∑w∈Wλwz~jh+1,wS tep 5.2 (set of freight terminal selected sequence with lack of empty vehicle). For the freight terminal where the number of load waiting to be arranged q~jh+1 is nonpositive, it means that there are empty vehicles enough not to continue to call the empty vehicles. Following the abovementioned principles, the freight terminals with lack of empty vehicle are ranked in descending order of the revenues ∑k∈Nαjkq~jkh+1. The set of freight terminal selected sequence with lack of empty vehicle J~h+1 is formed. We denote σ as the serial number of freight terminal with lack of empty vehicle in this set J~h+1 and σ^ as the total quantity of freight terminal with lack of empty vehicle. So the new set J~h+1 can be expressed in the form by(21)J~h+1=jσ∣jσ=j,j∈N,σ=1⋯,σ^S tep 5.3 (set of empty vehicles shortage amount selected sequence). We count the shortage amount of freight terminal with lack of empty vehicle in the set J~h+1. Further, the set of empty vehicles shortage amount selected sequence Q~jh+1 is obtained and is denoted as Q~jh+1=q~jh+1,j∈J~.S tep 5.4 (set of freight terminal selected sequence for dispatching empty vehicle). The destination terminaljis successively taken out from the set of freight terminal selected sequence with lack of empty vehicle J~h+1. Then, according to the vehicle routing matrix πw, the origin terminal i that has a link with destination terminal j is taken out. And, ranking the origin terminal i in ascending order of the distance, the set of freight terminal selected sequence for dispatching empty vehicle I~jh is formed. We denote θ as the serial number of freight terminal for dispatching empty vehicle in this set I~jh and σ^ as the total quantity of freight terminal for dispatching empty vehicle. So the new set I~jh can be expressed in the form by(22)I~jh=iθ∣πiθjw=1,iθ∈N,w∈W,j∈J~h+1,θ=1,⋯,θ^jhS tep 5.5 (empty vehicle reposition scheme). We choose the destination terminaljfrom the set of freight terminal selected sequence with lack of empty vehicle J~h+1 and the origin terminal i(θ) from the corresponding set of freight terminal selected sequence for dispatching empty vehicle I~jh. In accordance with vehicle distribution Z-i(θ)h at freight terminal for dispatching empty vehicle i(θ) and loads number q~jh+1 at freight terminal with lack of empty vehicle, we make use of the principle of priority arranging large capacity vehicle types to conduct the empty vehicle reposition as the following formula:(23)yjh,δθ=min⁡q~jh+1λδ,z-iθhδ,i∈I~jh(1) If there areq~jh+1/λδ≥z-i(θ)h(δ), then yjh,δ(θ)=z-i(θ)h(δ). Additionally, let q~jh+1=q~jh+1-λδz-i(θ)h(δ) and z-i(θ)h(δ)=0. If there are still some other vehicle types unarranged at freight terminal for dispatching empty vehicle i(θ), they are furtherly made by using equation (23). Otherwise, succeeding freight terminal i(θ) should be taken out from the set of freight terminal selected sequence for dispatching empty vehicle I~jh.(2) If there areq~jh+1/λδ≤z-i(θ)h(δ), then yjh,δ(θ)=q~jh+1/λδ. Besides, let q~jh+1=0 and z-i(θ)h(δ)=z-i(θ)h(δ)-q~jh+1/λδ. In addition, next destination stationj is taken out from set of freight terminal selected sequence with lack of empty vehicle J~h+1.Stage 6 (updating the vehicle distributions).
We use∑i∈Nxijh+1-tij,w+yijh+1-tij,wπijw=zjh+1,w to count the vehicle distributions Zjh+1,w in h period and return to Stage1.The approach of vehicle allocation with greedy heuristic procedure (VA-GHP) is shown in Figure4.Figure 4
Loaded vehicles dispatching and empty vehicle reposition scheme with greedy heuristic procedure.


## 6. Updating Solution Procedure by Improving Vehicle Allocation

In this section, we develop a method to update the initial solution generated by VA-GHP approach mentioned above. We name this method as updating solution procedure by improving vehicle allocation at local terminal. This method is abbreviated as USP-IVA in order to facilitate problem statement. Specially, the USP-IVA is composed of two subprocedures. The first is called updating solution by improving loaded vehicle dispatching scheme at single terminal and is abbreviated as US-ILVD. The second is denoted as updating solution by improving empty vehicle reposition scheme at single terminal and is abbreviated as US-IEVR. An overview of the framework of US-ILVD and US-IEVR is explained in Sections6.1 and 6.2.


### 6.1. Updating Solution by Improving Loaded Vehicle Dispatching Scheme at Single Terminal

It is not suitable to generate vehicle allocation scheme by the principle of priority arrangement for large load capacity vehicle types with VA-GHP approach for all freight terminals. In consideration of these disadvantages, we develop the solution improvement procedure called US-ILVD based on the initial solution. The concrete procedure of US-ILVD is as follows.Stage 1 (selecting the time period and freight terminal).
The time period and freight terminal are selected randomly from the initial vehicle allocation scheme generated by VA-GHP and denoted asi°,h°∣i°∈N,h°∈H.Stage 2 (updating loaded vehicle dispatching scheme of specific freight terminal and time period).
⁡
S tep 2.1 (adjusting the set of vehicle types selection sequence). According to the vehicle distribution zi°h°,w in terminal i° at h° time period, the new set of vehicle types selection sequence W⌢i°h° is obtained with ranking the vehicle types in random order. We denote δ as serial number of the vehicle type in this new set and δ^ as the total quantity of available vehicle types. Thus, W⌢i°h° can be recorded as(24)W⌢i°h°=w⌢i°h°δ∣w⌢i°h°δ=w,w∈W,δ=1,⋯,δ^S tep 2.2 (adjusting loaded vehicles dispatching scheme of transport routing). The transport routing si°h°(ξ) is successively chosen from the set of transport routing selection sequence Si°h° in order of the revenue from high to low. According to the vehicle routing matrix πw, the available vehicle types w⌢i°h°(δ) allowed to travel on the transport routing si°h°(ξ) is chosen from the new vehicle type selection sequence W⌢i°h°. And the fleet size z⌢i°h°(δ) of corresponding vehicle types is obtained from fleet size selection sequence Z⌢i°h°. We use the formula (25) to ensure loaded vehicles dispatching scheme:(25)xi°h°,δξ=min⁡qi°h°ξλδ,z⌢i°h°δS tep 2.3 (adjusting available loads number and available fleet size for transport routing). (1) If there are qi°h°(ξ)/λδ≥z⌢i°h°(δ), then xi°h°,δ(ξ)=z⌢i°h°(δ). Besides, let qi°h°(ξ)=qi°h°(ξ)-λδz⌢i°h°(δ) and z⌢i°h°(δ)=0; turn to Step 2.4.(2) If there areqi°h°(ξ)/λδ≤z⌢i°h°(δ), then xi°h°,δ(ξ)=qi°h°(ξ)/λδ. Besides, let qi°h°(ξ)=0 and z⌢i°h°(δ)=z⌢i°h°(δ)-qi°h°(ξ)/λδ; turn to Step2.5.S tep 2.4 (adjusting three selection sequence sets). We recount the revenue αijqi°h°(ξ) for each transport routing. The new set of transport routing selection sequence S-i°h° is updated with ranking the transport routing in descending order of the revenue. Then after removing the vehicle types w⌢i°h°(δ) used up from vehicle selection sequence set W⌢i°h°, the new set of vehicle type selection sequence W⌢⌢i°h° can be given. Further removing the fleet size z⌢i°h°(δ) of the vehicle types w⌢i°h°(δ) from fleet size selection sequence set Z⌢i°h°, the new set of fleet size selection sequence Z⌢⌢i°h° can be obtained.S tep 2.5(adjusting two selection sequence sets). The new set of transport routing selection sequence S-i°h° is given after removing the transport routing si°h°(ξ) from Si°h°. The new set of fleet sizing selection sequence Z⌢⌢i°h° is obtained after updating the corresponding fleet size z⌢i°h°(δ) of vehicle type w⌢i°h°(δ) in the set of fleet sizing selection sequence Z⌢i°h°.Stage 3 (determining loaded vehicle dispatching scheme of other freight terminals).
We work out the loaded vehicle dispatching scheme at other freight terminals in theh°th time period by the VA-GHP approach. And then the vehicle distributions at each freight terminal in the (h°+1)th time period are updated.Stage 4 (determining loaded vehicle dispatching and empty vehicle reposition scheme for remaining time periods).
The loaded vehicle dispatching and empty vehicle reposition scheme is obtained by VA-GHP approach for remaining time periods after theh°th time period. The updating solution by improving loaded vehicle dispatching scheme at single terminal is shown in Figure 5.Figure 5
Updating solutions by improving loaded vehicle dispatching scheme at single freight terminal.


### 6.2. Updating Solution by Improving Empty Vehicle Reposition Scheme at Single Terminal

On the basis of the improved solution obtained by US-ILVD, secondary improvement procedure called US-IEVR is further carried out. The concrete procedure of US-IEVR is as follows.Stage 1 (selecting the time period and freight terminal).
The time periodh∗ is selected randomly from planning horizon. And then the destination freight terminal j∗ is chosen randomly from freight terminal selected sequence with lack of empty vehicle J~h∗. They are denoted as j∗,h∗∣j∗∈N,h∗∈H.Stage 2 (adjusting the set of freight terminal selected sequence for dispatching empty vehicle).
According to the vehicle routing matrixπw, the origin terminal i∗ that has a link with destination terminal j∗ is taken out. And, ranking the origin terminal i∗ in ascending order of the distance, the set of freight terminal selected sequence for dispatching empty vehicle I~j∗h∗-1 is obtained. We denote i∗(θ) as the serial number of freight terminal for dispatching empty vehicle in this set I~j∗h∗-1.Stage 3 (adjusting empty vehicles reposition scheme).
We choose successively the origin terminali∗(θ) from the corresponding set of freight terminal selected sequence for dispatching empty vehicle I~j∗h∗-1. In accordance with vehicle distribution Z⌢i∗(θ)h∗-1 at freight terminal for dispatching empty vehicle i∗(θ) and loads number q~j∗h∗ at freight terminal with lack of empty vehicle j∗, we make use of the principle of randomly arranging vehicle types to generate the empty vehicle reposition scheme as the following formula:(26)yj∗h-1,δθ=min⁡q~j∗h∗λδ,z⌢i∗θh∗-1δ,i∗∈I~j∗h∗-1
(1) If there areq~jh∗/λδ≥z⌢i∗(θ)h∗-1(δ), yj∗h-1,δ(θ)=z⌢i∗(θ)h∗-1(δ). Additionally, let q~j∗h∗=q~j∗h∗-λδz⌢i∗(θ)h∗-1(δ) and z⌢i∗(θ)h∗-1(δ)=0. If there are still some other vehicle types in freight terminal for dispatching empty vehicle i∗(θ), we generate further the empty vehicle reposition scheme by equation (26). Otherwise, the next freight terminal i∗(θ) is taken out from the set of freight terminal selected sequence for dispatching empty vehicle I~j∗h∗-1.
(2) If there areq~jh∗/λδ≤z⌢i∗(θ)h∗-1(δ),yj∗h∗-1,δ(θ)=q~j∗h∗/λδ. Besides, let q~j∗h∗=0 and z⌢i∗(θ)h∗-1(δ)=z⌢i∗(θ)h∗-1(δ)-q~j∗h∗/λδ. And the procedures are stopped.Stage 4 (determining empty vehicle reposition scheme of other freight terminals).
We work out the empty vehicle reposition scheme at other freight terminals in theh∗th time period by the VA-GHP approach. And then the vehicle distributions at each freight terminal in the (h∗+1)th time period are updated.Stage 5 (determining loaded vehicle dispatching and empty vehicle reposition scheme for remaining time periods).
The loaded vehicle dispatching and empty vehicle reposition scheme is obtained by VA-GHP approach for remaining time periods after theh∗th time period. The updating solution by improving empty vehicle reposition scheme at single terminal is shown in Figure 6.Figure 6
Updating solutions by improving empty vehicle reposition scheme at single freight terminal.


## 6.1. Updating Solution by Improving Loaded Vehicle Dispatching Scheme at Single Terminal

It is not suitable to generate vehicle allocation scheme by the principle of priority arrangement for large load capacity vehicle types with VA-GHP approach for all freight terminals. In consideration of these disadvantages, we develop the solution improvement procedure called US-ILVD based on the initial solution. The concrete procedure of US-ILVD is as follows.Stage 1 (selecting the time period and freight terminal).
The time period and freight terminal are selected randomly from the initial vehicle allocation scheme generated by VA-GHP and denoted asi°,h°∣i°∈N,h°∈H.Stage 2 (updating loaded vehicle dispatching scheme of specific freight terminal and time period).
⁡
S tep 2.1 (adjusting the set of vehicle types selection sequence). According to the vehicle distribution zi°h°,w in terminal i° at h° time period, the new set of vehicle types selection sequence W⌢i°h° is obtained with ranking the vehicle types in random order. We denote δ as serial number of the vehicle type in this new set and δ^ as the total quantity of available vehicle types. Thus, W⌢i°h° can be recorded as(24)W⌢i°h°=w⌢i°h°δ∣w⌢i°h°δ=w,w∈W,δ=1,⋯,δ^S tep 2.2 (adjusting loaded vehicles dispatching scheme of transport routing). The transport routing si°h°(ξ) is successively chosen from the set of transport routing selection sequence Si°h° in order of the revenue from high to low. According to the vehicle routing matrix πw, the available vehicle types w⌢i°h°(δ) allowed to travel on the transport routing si°h°(ξ) is chosen from the new vehicle type selection sequence W⌢i°h°. And the fleet size z⌢i°h°(δ) of corresponding vehicle types is obtained from fleet size selection sequence Z⌢i°h°. We use the formula (25) to ensure loaded vehicles dispatching scheme:(25)xi°h°,δξ=min⁡qi°h°ξλδ,z⌢i°h°δS tep 2.3 (adjusting available loads number and available fleet size for transport routing). (1) If there are qi°h°(ξ)/λδ≥z⌢i°h°(δ), then xi°h°,δ(ξ)=z⌢i°h°(δ). Besides, let qi°h°(ξ)=qi°h°(ξ)-λδz⌢i°h°(δ) and z⌢i°h°(δ)=0; turn to Step 2.4.(2) If there areqi°h°(ξ)/λδ≤z⌢i°h°(δ), then xi°h°,δ(ξ)=qi°h°(ξ)/λδ. Besides, let qi°h°(ξ)=0 and z⌢i°h°(δ)=z⌢i°h°(δ)-qi°h°(ξ)/λδ; turn to Step2.5.S tep 2.4 (adjusting three selection sequence sets). We recount the revenue αijqi°h°(ξ) for each transport routing. The new set of transport routing selection sequence S-i°h° is updated with ranking the transport routing in descending order of the revenue. Then after removing the vehicle types w⌢i°h°(δ) used up from vehicle selection sequence set W⌢i°h°, the new set of vehicle type selection sequence W⌢⌢i°h° can be given. Further removing the fleet size z⌢i°h°(δ) of the vehicle types w⌢i°h°(δ) from fleet size selection sequence set Z⌢i°h°, the new set of fleet size selection sequence Z⌢⌢i°h° can be obtained.S tep 2.5(adjusting two selection sequence sets). The new set of transport routing selection sequence S-i°h° is given after removing the transport routing si°h°(ξ) from Si°h°. The new set of fleet sizing selection sequence Z⌢⌢i°h° is obtained after updating the corresponding fleet size z⌢i°h°(δ) of vehicle type w⌢i°h°(δ) in the set of fleet sizing selection sequence Z⌢i°h°.Stage 3 (determining loaded vehicle dispatching scheme of other freight terminals).
We work out the loaded vehicle dispatching scheme at other freight terminals in theh°th time period by the VA-GHP approach. And then the vehicle distributions at each freight terminal in the (h°+1)th time period are updated.Stage 4 (determining loaded vehicle dispatching and empty vehicle reposition scheme for remaining time periods).
The loaded vehicle dispatching and empty vehicle reposition scheme is obtained by VA-GHP approach for remaining time periods after theh°th time period. The updating solution by improving loaded vehicle dispatching scheme at single terminal is shown in Figure 5.Figure 5
Updating solutions by improving loaded vehicle dispatching scheme at single freight terminal.


## 6.2. Updating Solution by Improving Empty Vehicle Reposition Scheme at Single Terminal

On the basis of the improved solution obtained by US-ILVD, secondary improvement procedure called US-IEVR is further carried out. The concrete procedure of US-IEVR is as follows.Stage 1 (selecting the time period and freight terminal).
The time periodh∗ is selected randomly from planning horizon. And then the destination freight terminal j∗ is chosen randomly from freight terminal selected sequence with lack of empty vehicle J~h∗. They are denoted as j∗,h∗∣j∗∈N,h∗∈H.Stage 2 (adjusting the set of freight terminal selected sequence for dispatching empty vehicle).
According to the vehicle routing matrixπw, the origin terminal i∗ that has a link with destination terminal j∗ is taken out. And, ranking the origin terminal i∗ in ascending order of the distance, the set of freight terminal selected sequence for dispatching empty vehicle I~j∗h∗-1 is obtained. We denote i∗(θ) as the serial number of freight terminal for dispatching empty vehicle in this set I~j∗h∗-1.Stage 3 (adjusting empty vehicles reposition scheme).
We choose successively the origin terminali∗(θ) from the corresponding set of freight terminal selected sequence for dispatching empty vehicle I~j∗h∗-1. In accordance with vehicle distribution Z⌢i∗(θ)h∗-1 at freight terminal for dispatching empty vehicle i∗(θ) and loads number q~j∗h∗ at freight terminal with lack of empty vehicle j∗, we make use of the principle of randomly arranging vehicle types to generate the empty vehicle reposition scheme as the following formula:(26)yj∗h-1,δθ=min⁡q~j∗h∗λδ,z⌢i∗θh∗-1δ,i∗∈I~j∗h∗-1
(1) If there areq~jh∗/λδ≥z⌢i∗(θ)h∗-1(δ), yj∗h-1,δ(θ)=z⌢i∗(θ)h∗-1(δ). Additionally, let q~j∗h∗=q~j∗h∗-λδz⌢i∗(θ)h∗-1(δ) and z⌢i∗(θ)h∗-1(δ)=0. If there are still some other vehicle types in freight terminal for dispatching empty vehicle i∗(θ), we generate further the empty vehicle reposition scheme by equation (26). Otherwise, the next freight terminal i∗(θ) is taken out from the set of freight terminal selected sequence for dispatching empty vehicle I~j∗h∗-1.
(2) If there areq~jh∗/λδ≤z⌢i∗(θ)h∗-1(δ),yj∗h∗-1,δ(θ)=q~j∗h∗/λδ. Besides, let q~j∗h∗=0 and z⌢i∗(θ)h∗-1(δ)=z⌢i∗(θ)h∗-1(δ)-q~j∗h∗/λδ. And the procedures are stopped.Stage 4 (determining empty vehicle reposition scheme of other freight terminals).
We work out the empty vehicle reposition scheme at other freight terminals in theh∗th time period by the VA-GHP approach. And then the vehicle distributions at each freight terminal in the (h∗+1)th time period are updated.Stage 5 (determining loaded vehicle dispatching and empty vehicle reposition scheme for remaining time periods).
The loaded vehicle dispatching and empty vehicle reposition scheme is obtained by VA-GHP approach for remaining time periods after theh∗th time period. The updating solution by improving empty vehicle reposition scheme at single terminal is shown in Figure 6.Figure 6
Updating solutions by improving empty vehicle reposition scheme at single freight terminal.


## 7. Hybrid Simulated Annealing Algorithm Combining VA-GHP&USP-IVA

In this section, we develop a hybrid simulated annealing heuristic to solve MHFS-FSD using the VA-GHP approach and USP-IVA method described previously in Sections5 and 6. In this proposed algorithm, we use the VA-GHP procedure to generate initial solutions and the USP-IVA procedure to improve the solutions. We incorporate VA-GHP and USP-IVA into simulated annealing algorithm. So the hybrid simulated annealing heuristic is proposed. The detailed steps of hybrid simulated annealing heuristics (HSAH) are as follows.


### 7.1. Acceptance Measure of Solutions

(  1) Acceptance Measure for US-ILVD. Let ψ∗=Xh,Yh,Zh∣h=1,⋯,K denote the solution accepted currently and ψUS-ILVD(i°,h°)=ψ∗→X⌢i°h° express the new solution obtained by US-ILVD of single terminal (i°,h°). We let FUS-ILVD+i°,h° be the incremental cost as follows:(27)FUS-ILVD+i°,h°=FψUS-ILVDi°,h°-Fψ∗The acceptance measure of new solutionψUS-ILVD(i°,h°) is written as follows:(28)PψUS-ILVDi°,h°=1,ifFUS-ILVD+i°,h°≤0.exp⁡-FUS-ILVD+i°,h°t,ifFUS-ILVD+i°,h°>0.whereP(ψUS-ILVD(i°,h°)) denotes the acceptance probability of new solution ψUS-ILVD(i°,h°).(  2) Acceptance Measure for US-IEVR. Let ψ∗=Xh,Yh,Zh∣h=1,⋯,K denote the solution accepted currently and ψUS-IEVR(j∗,h∗)=ψ∗→Y⌢j∗h∗ express the solution obtained by US-IEVR of single terminal (j∗,h∗). Let us define FUS-IEVR+(j∗,h∗) as the incremental cost. We have that(29)FUS-IEVR+j∗,h∗=FψUS-IEVRj∗,h∗-Fψ∗We can form the acceptance measure of new solutionψUS-IEVR(j∗,h∗) as follows:(30)PψUS-IEVRj∗,h∗=1,ifFUS-IEVR+j∗,h∗≤0.exp⁡-FUS-IEVR+j∗,h∗tifFUS-IEVR+j∗,h∗>0.whereP(ψUS-IEVR(j∗,h∗)) is the acceptance probability of new solution ψUS-IEVR(j∗,h∗).


### 7.2. Control Parameters of the Heuristics

⁡(  1) Initial Temperature Parameter t 0. In order to return a high quality final solution at the end of the algorithm iteration, we must give a large enough value to t0. But the exorbitant initial temperature can cause the sharp growth of CPU run time so that makes the algorithm lose its feasibility. So we let t0=2×F(ψ∗), where F(ψ∗) is the objective function value obtained by taking the vehicles allocation scheme ψ=X1,Y1,Z1,…,XK,YK,ZK into objective function of MHFS-FSD model.(  2) Termination Criterion. Here we give a very small positive number tf as the termination temperature of hybrid simulated annealing heuristics. During the iterative procedure of the algorithm, the temperature control parameter t is gradually reduced. When the temperature control parameter t is less than tf, the algorithm will be terminated.(  3) Attenuation Function of Temperature Control Parameter. In order to avoid overlong running times in each iteration, the decrement of temperature control parametertshould take a small value. Here we have the attenuation function of temperature control parameter as follows:(31)tk+1=α·tk,k=0,1,…,0<α<1wheretk and tk+1 are, respectively, the temperature control parameter in kth iteration and (k+1)th iteration, α is the attenuation coefficient, and we assume α to be 0.9.(  4) Solution Transformation Times. The acceptance probability of new solution will reduce with the decrease of the temperature control parameter tk. Consequently, solution transformation times Lk should gradually increase when temperature control parameter tk drops; that is, tk↓⟹Lk↑. Here we have the dynamic solution transformation times based on algorithm iterative procedure:(32)Lk+1=β·Lk,k=0,1,…,β>1whereLk and Lk+1 are, respectively, the solution transformation times in kth iteration and (k+1)th iteration and β is the increasing rate. We assume β to be 1.1 and L0 to be a number between 50~100.


### 7.3. Heuristic Procedure

For illustration purpose, the hybrid simulated annealing heuristic combining greedy heuristic procedure for initial solution and improving vehicle allocation for updating solution is abbreviated as HSAH-GHP&IMA. Figure7 depicts a flowchart of the hybrid simulated annealing heuristics. The heuristic procedure can be described as follows.Figure 7
Flow diagram of hybrid simulated annealing heuristics.Stage 1 (initial scheduling scheme generated by GHP).
Letψ=X1,Y1,Z1,…,XK,YK,ZK denote the vehicle allocation scheme of whole planning horizon obtained by VA-GHP. And then the objective function value F(ψ) can be obtained by taking it into the objective function of MHFS-FSD model. Besides, we let ψ∗=ψ.Stage 2 (setting control parameters).
We set the initial temperaturet0, termination temperature tf, attenuation coefficient of temperature control parameter α, initial solution transformation times L0, and increasing rate of solution transformation times β. And then let t=t0, L=L0; turn to Step2.Stage 3 (iteration based on cooling schedule).
⁡
S tep 1 (measuring termination condition). In case t>tf, the algorithm procedure turns to Step2. Otherwise, move to Step4.S tep 2(implementing US-IVA procedure). Let i=1.S tep 2.1(implementing US-ILVD subprocedure). After choosing randomly the adjustment freight terminal and time period (i°,h°), we use US-ILVD procedure to improve solution ψ∗ so that the new solution ψUS-ILVD(i°,h°) is obtained. The FUS-ILVD+(i°,h°) value is calculated by equation (27). If FUS-ILVD+(i°,h°)≤0, we accept the new solution and let ψ∗=ψUS-ILVD(i°,h°). If a random number R between 0 and 1 is generated and exp⁡(-FUS-ILVD+(i°,h°)/t)≥R, we accept the new solution and let ψ∗=ψUS-ILVD(i°,h°). Otherwise, we reject to accept the new solution and ψ∗ remain unchanged. The procedure goes to Step2.2.S tep 2.2(implementing US- IEVR subprocedure). After choosing randomly the adjustment freight terminal and time period (j∗,h∗), we use US-IEVR procedure to improve solution ψ∗ so that the new solution ψUS-IEVR(j∗,h∗) is obtained. The FUS-IEVR+(j∗,h∗) value is calculated by equation (29). If FUS-IEVR+(j∗,h∗)≤0, we accept the new solution and let ψ∗=ψUS-IEVR(j∗,h∗). If a random number R between 0 and 1 is generated and exp⁡(-FUS-IEVR+(j∗,h∗)/t)≥R, we accept the new solution and let ψ∗=ψUS-IEVR(j∗,h∗). Otherwise, we reject to accept the new solution and ψ∗ remain unchanged. The procedure goes to Step2.3.S tep 2.3(counting solution transformation times). Let i=i+1. If there is i≤L, turn to Step2.1; else go to Step2.4.S tep 2.4 (adjusting temperature and solution transformation times). Let t=t×α and L=L×β, and go to Step1.Stage 4 (obtaining final solution).
When the termination criterion is met, the final solutionψ∗ can be obtained.


## 7.1. Acceptance Measure of Solutions

(  1) Acceptance Measure for US-ILVD. Let ψ∗=Xh,Yh,Zh∣h=1,⋯,K denote the solution accepted currently and ψUS-ILVD(i°,h°)=ψ∗→X⌢i°h° express the new solution obtained by US-ILVD of single terminal (i°,h°). We let FUS-ILVD+i°,h° be the incremental cost as follows:(27)FUS-ILVD+i°,h°=FψUS-ILVDi°,h°-Fψ∗The acceptance measure of new solutionψUS-ILVD(i°,h°) is written as follows:(28)PψUS-ILVDi°,h°=1,ifFUS-ILVD+i°,h°≤0.exp⁡-FUS-ILVD+i°,h°t,ifFUS-ILVD+i°,h°>0.whereP(ψUS-ILVD(i°,h°)) denotes the acceptance probability of new solution ψUS-ILVD(i°,h°).(  2) Acceptance Measure for US-IEVR. Let ψ∗=Xh,Yh,Zh∣h=1,⋯,K denote the solution accepted currently and ψUS-IEVR(j∗,h∗)=ψ∗→Y⌢j∗h∗ express the solution obtained by US-IEVR of single terminal (j∗,h∗). Let us define FUS-IEVR+(j∗,h∗) as the incremental cost. We have that(29)FUS-IEVR+j∗,h∗=FψUS-IEVRj∗,h∗-Fψ∗We can form the acceptance measure of new solutionψUS-IEVR(j∗,h∗) as follows:(30)PψUS-IEVRj∗,h∗=1,ifFUS-IEVR+j∗,h∗≤0.exp⁡-FUS-IEVR+j∗,h∗tifFUS-IEVR+j∗,h∗>0.whereP(ψUS-IEVR(j∗,h∗)) is the acceptance probability of new solution ψUS-IEVR(j∗,h∗).


## 7.2. Control Parameters of the Heuristics

⁡(  1) Initial Temperature Parameter t 0. In order to return a high quality final solution at the end of the algorithm iteration, we must give a large enough value to t0. But the exorbitant initial temperature can cause the sharp growth of CPU run time so that makes the algorithm lose its feasibility. So we let t0=2×F(ψ∗), where F(ψ∗) is the objective function value obtained by taking the vehicles allocation scheme ψ=X1,Y1,Z1,…,XK,YK,ZK into objective function of MHFS-FSD model.(  2) Termination Criterion. Here we give a very small positive number tf as the termination temperature of hybrid simulated annealing heuristics. During the iterative procedure of the algorithm, the temperature control parameter t is gradually reduced. When the temperature control parameter t is less than tf, the algorithm will be terminated.(  3) Attenuation Function of Temperature Control Parameter. In order to avoid overlong running times in each iteration, the decrement of temperature control parametertshould take a small value. Here we have the attenuation function of temperature control parameter as follows:(31)tk+1=α·tk,k=0,1,…,0<α<1wheretk and tk+1 are, respectively, the temperature control parameter in kth iteration and (k+1)th iteration, α is the attenuation coefficient, and we assume α to be 0.9.(  4) Solution Transformation Times. The acceptance probability of new solution will reduce with the decrease of the temperature control parameter tk. Consequently, solution transformation times Lk should gradually increase when temperature control parameter tk drops; that is, tk↓⟹Lk↑. Here we have the dynamic solution transformation times based on algorithm iterative procedure:(32)Lk+1=β·Lk,k=0,1,…,β>1whereLk and Lk+1 are, respectively, the solution transformation times in kth iteration and (k+1)th iteration and β is the increasing rate. We assume β to be 1.1 and L0 to be a number between 50~100.


## 7.3. Heuristic Procedure

For illustration purpose, the hybrid simulated annealing heuristic combining greedy heuristic procedure for initial solution and improving vehicle allocation for updating solution is abbreviated as HSAH-GHP&IMA. Figure7 depicts a flowchart of the hybrid simulated annealing heuristics. The heuristic procedure can be described as follows.Figure 7
Flow diagram of hybrid simulated annealing heuristics.Stage 1 (initial scheduling scheme generated by GHP).
Letψ=X1,Y1,Z1,…,XK,YK,ZK denote the vehicle allocation scheme of whole planning horizon obtained by VA-GHP. And then the objective function value F(ψ) can be obtained by taking it into the objective function of MHFS-FSD model. Besides, we let ψ∗=ψ.Stage 2 (setting control parameters).
We set the initial temperaturet0, termination temperature tf, attenuation coefficient of temperature control parameter α, initial solution transformation times L0, and increasing rate of solution transformation times β. And then let t=t0, L=L0; turn to Step2.Stage 3 (iteration based on cooling schedule).
⁡
S tep 1 (measuring termination condition). In case t>tf, the algorithm procedure turns to Step2. Otherwise, move to Step4.S tep 2(implementing US-IVA procedure). Let i=1.S tep 2.1(implementing US-ILVD subprocedure). After choosing randomly the adjustment freight terminal and time period (i°,h°), we use US-ILVD procedure to improve solution ψ∗ so that the new solution ψUS-ILVD(i°,h°) is obtained. The FUS-ILVD+(i°,h°) value is calculated by equation (27). If FUS-ILVD+(i°,h°)≤0, we accept the new solution and let ψ∗=ψUS-ILVD(i°,h°). If a random number R between 0 and 1 is generated and exp⁡(-FUS-ILVD+(i°,h°)/t)≥R, we accept the new solution and let ψ∗=ψUS-ILVD(i°,h°). Otherwise, we reject to accept the new solution and ψ∗ remain unchanged. The procedure goes to Step2.2.S tep 2.2(implementing US- IEVR subprocedure). After choosing randomly the adjustment freight terminal and time period (j∗,h∗), we use US-IEVR procedure to improve solution ψ∗ so that the new solution ψUS-IEVR(j∗,h∗) is obtained. The FUS-IEVR+(j∗,h∗) value is calculated by equation (29). If FUS-IEVR+(j∗,h∗)≤0, we accept the new solution and let ψ∗=ψUS-IEVR(j∗,h∗). If a random number R between 0 and 1 is generated and exp⁡(-FUS-IEVR+(j∗,h∗)/t)≥R, we accept the new solution and let ψ∗=ψUS-IEVR(j∗,h∗). Otherwise, we reject to accept the new solution and ψ∗ remain unchanged. The procedure goes to Step2.3.S tep 2.3(counting solution transformation times). Let i=i+1. If there is i≤L, turn to Step2.1; else go to Step2.4.S tep 2.4 (adjusting temperature and solution transformation times). Let t=t×α and L=L×β, and go to Step1.Stage 4 (obtaining final solution).
When the termination criterion is met, the final solutionψ∗ can be obtained.


## 8. Numerical Experiments

In this section, we try to evaluate the quality of the hybrid simulated annealing algorithm combining VA-GHP&USP-IVA proposed in previous section to solve the MHFS-FSD in terms of traditional measure such as objective function and execution time. The scenarios for MHFS-FSD are described in Section8.1 and three sets of experiments are devised to test the effectiveness of the hybrid simulated annealing algorithm proposed. Section 8.2 reports the numerical results.


### 8.1. Scenario Settings of the Experiments

To design the experiments scheme, the scenarios settings need to be firstly set up. This section describes the data used in the numerical testing of the models. We generated the dataset based on the freight enterprise in Henan Province, China. However, the data have some inconsistencies caused by multiple types of the vehicles, differences in financial accounting system of various corporations, and so on. Therefore the data has to be cleaned.Here we chose a certain range from 40 RMB to 60 RMB per load as the revenue of loaded movement between various freight terminals. The loading capacity of various type vehicles is allowed the certain range from 20 units to 30 units per vehicle. The cost of empty reposition is allowed the certain range from 200 RMB to 300 RMB per vehicle. We use the cost range from 2000 RMB to 3000 RMB per vehicles as purchase or lease cost of various vehicles.To evaluate the performance of the heuristics proposed, we generate randomly three sets of problem instances with different size of freight terminal and time periods, considering different number of vehicles types. Three sets of experiments are oriented to compare the performances of the hybrid simulated annealing algorithm with the CPLEX solver and VA-GHP approach.In the first set of contrast experiment, there are the heterogeneous fleet with three types of vehicles and five time periods. The number of the freight terminal is from 5 to 10. In the second set of contrast experiment, there are the heterogeneous fleet with five types of vehicles and ten time periods. The number of the freight terminal is from 11 to 20. In the third set of contrast experiment, there are the heterogeneous fleet with ten types of vehicles and thirty time periods. The number of the freight terminal is from 21 to 50.The test work has been done using a computer with Intel i5-5200CPU (2.2 GHz). The hybrid simulated annealing algorithm combining VA-GHP&USP-IVA has been coded in Microsoft Visual C++.


### 8.2. Performance Evaluation

In this section the computational experiments are performed using three methods to solve the MHFS-FSD model: HSAH program, GHP approach, and CPLEX solver. And we use two measures to evaluate the performance. The first one is the OPT, which is the value of the objective function obtained by the CPLEX solver, GHP approach, and HSAH program. OPT as one of the major criteria in assessing the performance of the HSAH, GHP, and CPLEX represents the profit generated by revenues for dispatching loaded vehicles, repositioning empty vehicles, and ownership costs for owning vehicle in planning horizon. The second measure of performance is the CUP time to run CPLEX solver, the GHP approach, and HSAH program for MHFS-FSD model. We aim at quantifying the benefits of the proposed method in two aspects: (1) the efficiency of the proposed method executed with no time limit imposed and (2) the efficiency of the proposed method executed with time limit imposed.(  1) Measure of Performance with No CPU Time Limit Imposed. In this section, we focus on the numerical studies on the performance of GHP approach, HSAH program and CPLEX solver executed with no CPU time limit imposed. We compare the performance of these three methods in two aspects: solution quality and computational efficiency. In Tables 1–3 we present the results obtained by solving randomly generated MHFS-FSD problems. We compare the solutions obtained by solving the MHFS-FSD problems using GHP approach, HSAH program, and CPLEX solver.Table 1
Performance for CPLEX, GHP, and HSAH applied to small size problem with no time limit imposed.
Problem size (freight terminal × vehicle types × time period) CPLEX GHP HSAH OS CPU time OFV CPU time RAT GAP OFV CPU time RAT GAP (RMB) (second) (RMB) (second) (%) (second) (RMB) (second) (%) (second) 5 × 3 × 5 219586 214.5 165894 13.6 24.45 200.9 199473 33.6 9.16 180.9 6 × 3 × 5 225465 235.7 184759 15.7 18.05 220.0 209484 37.5 7.09 198.2 7 × 3 × 5 230187 356.2 163028 13.9 29.18 342.3 190583 36.6 17.21 319.6 8 × 3 × 5 246586 248.3 190324 21.6 22.82 226.7 224756 45.6 8.85 202.7 9 × 3 × 5 250829 321.2 186206 25.8 25.76 295.4 208573 50.4 16.85 270.8 10 × 3 × 5 276585 314.5 201478 23.1 27.16 291.4 235859 55.6 14.72 258.9Table 2
Performance for CPLEX, GHP, and HSAH applied to medium size problem with no time limit imposed.
Problem size (freight terminal × vehicle types × time period) CPLEX GHP HSAH OS CPU time OFV CPU time RAT GAP OFV CPU time RAT GAP (RMB) (second) (RMB) (second) (%) (second) (RMB) (second) (%) (second) 11 × 5 × 10 727865 389.5 512034 89.7 29.65 299.8 708543 196.7 2.65 192.8 12 × 5 × 10 767289 397.2 536401 91.4 30.09 305.8 741478 183.4 3.36 213.8 13 × 5 × 10 876537 405.7 641963 103.4 26.76 302.3 835317 209.9 4.70 195.8 14 × 5 × 10 782685 412.3 541069 113.5 30.87 298.8 729294 220.3 6.82 192.0 15 × 5 × 10 813858 434.2 571648 73.8 29.76 360.4 787625 177.8 3.22 256.4 16 × 5 × 10 882385 451.8 642917 112.5 27.13 339.3 837674 229.6 5.07 222.2 17 × 5 × 10 926048 446.9 702541 100.9 24.14 346.0 892573 214.2 3.61 232.7 18 × 5 × 10 949874 467.1 731529 85.2 22.99 381.9 932074 192.1 1.87 275.0 19 × 5 × 10 998375 482.6 768359 78.4 23.04 404.2 948927 187.4 4.95 295.2 20 × 5 × 10 1010948 513.7 763542 104.3 24.47 409.4 978375 210.6 3.22 303.1Table 3
Performance for CPLEX, GHP, and HSAH applied to large size problem with no time limit imposed.
Problem size (freight terminal × vehicle types × time period) CPLEX GHP HSAH OS CPU time OFV CPU time RAT GAP OFV CPU time RAT GAP (RMB) (second) (RMB) (second) (%) (second) (RMB) (second) (%) (second) 21 × 10 × 30 3641457 1046.6 3162419 351.7 13.16 694.9 3592236 631.5 1.35 415.1 22 × 10 × 30 3960247 1153.7 3375864 363.2 14.76 790.5 3890437 654.2 1.76 499.5 23 × 10 × 30 4093746 1241.3 3392618 389.1 17.13 852.2 3962476 691.0 3.21 550.3 24 × 10 × 30 4172590 1306.9 3521607 412.5 15.60 894.4 4063582 719.6 2.61 587.3 25 × 10 × 30 4395799 1388.1 3862453 434.7 12.13 953.4 4298741 767.9 2.21 620.2 26 × 10 × 30 4673825 1449.7 3971641 443.8 15.02 1005.9 4495476 779.0 3.82 670.7 27 × 10 × 30 4691793 1533.8 4071358 503.4 13.22 1030.4 4586179 839.1 2.25 694.7 28 × 10 × 30 4762133 1592.4 4193062 536.7 11.95 1055.7 4703167 864.2 1.24 728.2 29 × 10 × 30 4792619 1640.3 4083017 576.3 14.81 1064.0 4698731 899.6 1.96 740.7 30 × 10 × 30 4821617 1731.6 4262843 598.6 11.59 1133.0 4803164 904.3 0.38 827.3 31 × 10 × 30 4865329 1796.1 4228019 602.6 13.10 1193.5 4829316 932.5 0.74 863.6 32 × 10 × 30 4887691 1874.2 4353762 647.9 10.92 1226.3 4796317 966.7 1.87 907.5 33 × 10 × 30 4963720 1924.0 4464907 712.4 10.05 1211.6 4850691 1043.2 2.28 880.8 34 × 10 × 30 5120379 1989.6 4526143 726.8 11.61 1262.8 5006479 1132.5 2.22 857.1 35 × 10 × 30 5436826 2046.5 4593011 779.3 15.52 1267.2 5267984 1164.6 3.11 881.9 36 × 10 × 30 5761844 2165.3 4612938 852.4 19.94 1312.9 5706914 1254.2 0.95 911.1 37 × 10 × 30 5839764 2221.2 4685714 837.0 19.76 1384.2 5784639 1275.6 0.94 945.6 38 × 10 × 30 6032875 2302.2 4703519 954.3 22.04 1347.9 5859217 1324.8 2.88 977.4 39 × 10 × 30 6361448 2399.5 4794808 989.4 24.63 1410.1 6263941 1396.1 1.53 1003.4 40 × 10 × 30 6673280 2556.4 4862714 1032.6 27.13 1523.8 6612476 1443.2 0.91 1113.2 41 × 10 × 30 6805218 2826.5 5182913 1058.1 23.84 1768.4 6569031 1504.1 3.47 1322.4 42 × 10 × 30 6961760 2973.3 5403917 1129.8 22.38 1843.5 6795213 1575.4 2.39 1397.9 43 × 10 × 30 7134752 3205.6 5632051 1186.7 21.06 2018.9 6948217 1614.5 2.61 1591.1 44 × 10 × 30 7446389 3360.5 5789830 1056.3 22.25 2304.2 7164206 1653.2 3.79 1707.3 45 × 10 × 30 7623796 3586.6 5913517 1219.4 22.43 2367.2 7501286 1741.8 1.61 1844.8 46 × 10 × 30 7853219 3888.2 6242683 1234.7 20.51 2653.5 7632689 1813.6 2.81 2074.6 47 × 10 × 30 7905256 4169.6 6590263 1196.4 16.63 2973.2 7693618 1856.8 2.68 2312.8 48 × 10 × 30 8134962 4275.3 6832657 1324.7 16.01 2950.6 8047603 1943.6 1.07 2331.7 49 × 10 × 30 8416732 4453.6 7207149 1717.2 14.37 2736.4 8267798 2066.2 1.77 2387.4 50 × 10 × 30 8601795 4679.1 7421678 1536.8 13.72 3142.3 8396419 2164.7 2.39 2514.4CPLEX solver is executed with no time limit imposed. The parameters used for the HSAH are as follows. The initial temperature is set byt0=2×F(ψ∗), where F(ψ∗) is the objective function value obtained by taking the initial solution solved with VA-IGP into objective function of MHFS-FSD model. The termination temperature is set to 1. The attenuation coefficient of temperature control parameter is assumed to be 0.9. The initial solution transformation frequency is allowed 80. The increasing rate of the solution transformation frequency is assumed to be 1.1.In order to facilitate problem analysis, we introduce two analytical indicators. The first denoted as RAT represents the deviation rate from the optimal solution. Let OFV be the objective function value with GHP approach and HSAH program. And let OS be optimal solution obtained with CPLEX solver. The RAT is calculated by(33)RAT=OS-OFVOSThe second is the CPU time gap denoted as GAP. It is calculated by the difference between CUP time of running the HSAH program or GHP approach and that of CPLEX solver. The CPU times of running HSAH program, GHP approach, and CPLEX solver are, respectively, denoted as CPUT.HSAH, CPU.GHP, and CPUT.CPLEX. The GAP is calculated by(34)GAP=CPUT.CPLEX-CPUT.HSAHGAP=CPUT.CPLEX-CPUT.GHPThe performance of small size problem is shown in Table1. Here we have 6 test problems. Each problem is, respectively, solved using the proposed GHP approach, HSAH program, and CPLEX solver. For the small size problem, the performance of GHP and HSAH is not very prominent when compared to CPLEX. And the performance of HSAH is better than GHP. The deviations of GHP from the optimal solution are the range from 18% to 30% while the deviations of HSAH from the optimal solution are the range from 7% to 20%. The CPU times to obtain the solutions are the range from 33.6 s to 55.6 s for HSAH, the range from 13.6 s to 23.1 s for GHP, and the range from 214.5 s to 356.2 s for CPLEX. The superiority of GHP and HSAH in computing time can be shown distinctly. The results of comparison are displayed in Figure 8.Figure 8
Results of comparison for small size problem.
(a)
Objective function value (b)
Deviation rate (c)
CPU running time (d)
CPU time gapThe performance of medium size problem is shown in Table2. Here we have 10 test problems. Each problem is respectively solved using the proposed GHP, HSAH, and CPLEX. For the medium size problem, the GHP and HSAH perform well when compared to CPLEX. The performance of HSAH is prominent than GHP. The deviations from the optimal solution are reduced to a reasonable range from 2% to 7% for HSAH and the deviations from the optimal solution range from 22% to 31% for GHP. One obvious observation is that the increase of CUP time to obtain the solutions using CPLEX begins to appear a trend of acceleration. But the rising step of CPU times to obtain the solutions using GHP and HSAH increases is relatively flat. So the superiority of GHP and HSAH in computation time is very obvious. The performance comparisons are shown in Figure 9.Figure 9
Results of comparison for medium size problem.
(a)
Objective function value (b)
Deviation rate (c)
CPU running time (d)
CPU time gapThe performance of large size problem is shown in Table3. Here we have 30 test problems. Each problem is respectively solved using the proposed GHP, HSAH, and CPLEX. For the large size problem, performance of solving the problem using the GHP and HSAH is more excellent than that of CPLEX. The deviations from the optimal solution for HSAH are limited to a very small range within 4% while the deviations from the optimal solution for the GHP are limited to range within 30%. But the CPU times for obtaining the solutions using CPLEX solver increase in a sharp speed due to the high memory requirements. The superiority of GHP and HSAH in computation time is particularly obvious. The performance comparisons are shown in Figure 10.Figure 10
Results of comparison for large size problem.
(a)
Objective function value (b)
Deviation rate (c)
CPU running time (d)
CPU time gap(  2) Measure of Performance with CPU Time Limit Imposed. In the experiment, we test the performance of the HSAH program and CPLEX solver with CPU time limit imposed. Since the performance of GHP approach is obviously inferior to that of HSAS program, GHP approach is no longer tested in the experiment. Here the HSAS program is fully executed. And let CPLEX solver be executed for the same CPU time. Then we compare the objective function value obtained by solving the MHFS-FSD problems using HSAH program and CPLEX solver. The tests are performed mainly in Tables 4–6.Table 4
Results for HSAH and CPLEX applied to small size problem with time limit imposed.
Problem size (freight terminal × vehicle types × time period) OFV (RMB) CPU time (second) CPLEX HSAH 5 × 3 × 5 53817 199473 33.6 6 × 3 × 5 57429 209484 37.5 7 × 3 × 5 61437 190583 36.6 8 × 3 × 5 68413 224756 45.6 9 × 3 × 5 73076 208573 50.4 10 × 3 × 5 78462 235859 55.6Table 5
Results for HSAH and CPLEX applied to medium size problem with time limit imposed.
Problem size (freight terminal × vehicle types × time period) OFV (RMB) CPU time (second) CPLEX HSAH 11 × 5 × 10 493276 708543 196.7 12 × 5 × 10 497625 741478 183.4 13 × 5 × 10 513410 835317 209.9 14 × 5 × 10 526683 729294 220.3 15 × 5 × 10 525914 787625 177.8 16 × 5 × 10 543728 837674 229.6 17 × 5 × 10 562807 892573 214.2 18 × 5 × 10 567436 932074 192.1 19 × 5 × 10 573614 948927 187.4 20 × 5 × 10 592460 978375 210.6Table 6
Results for HSAH and CPLEX applied to large size problem with time limit imposed.
Problem size (freight terminal × vehicle types × time period) OFV (RMB) CPU time (second) CPLEX HSAH 21 × 10 × 30 2561437 3592236 631.5 22 × 10 × 30 2746308 3890437 654.2 23 × 10 × 30 3012561 3962476 691.0 24 × 10 × 30 3348021 4063582 719.6 25 × 10 × 30 3585719 4298741 767.9 26 × 10 × 30 3602476 4495476 779.0 27 × 10 × 30 3680418 4586179 839.1 28 × 10 × 30 3764292 4703167 864.2 29 × 10 × 30 3791068 4698731 899.6 30 × 10 × 30 3822516 4803164 904.3 31 × 10 × 30 3849317 4829316 932.5 32 × 10 × 30 3886203 4796317 966.7 33 × 10 × 30 3965042 4850691 1043.2 34 × 10 × 30 4002837 5006479 1132.5 35 × 10 × 30 4125340 5267984 1164.6 36 × 10 × 30 4196205 5706914 1254.2 37 × 10 × 30 4273568 5784639 1275.6 38 × 10 × 30 4305217 5859217 1324.8 39 × 10 × 30 4362208 6263941 1396.1 40 × 10 × 30 4385426 6612476 1443.2 41 × 10 × 30 4412716 6569031 1504.1 42 × 10 × 30 4472930 6795213 1575.4 43 × 10 × 30 4502814 6948217 1614.5 44 × 10 × 30 4575356 7164206 1653.2 45 × 10 × 30 4598325 7501286 1741.8 46 × 10 × 30 4620819 7632689 1813.6 47 × 10 × 30 4673369 7693618 1856.8 48 × 10 × 30 4702923 8047603 1943.6 49 × 10 × 30 4785527 8267798 2066.2 50 × 10 × 30 4832509 8396419 2164.7From above three sets of contrast tests, the objective function values obtained by HSAH program are significantly better than those obtained by CPLEX solver with CPU time limit imposed.In sum, we can conclude that the HSAH is effective to obtain good solutions in relatively low run times for the MHFS-FSD. The results obtained indicate the benefits of using HSAH, helping MHFS-FSD to find better solutions in reasonable short CPU times, which are acceptable for practical applications.


## 8.1. Scenario Settings of the Experiments

To design the experiments scheme, the scenarios settings need to be firstly set up. This section describes the data used in the numerical testing of the models. We generated the dataset based on the freight enterprise in Henan Province, China. However, the data have some inconsistencies caused by multiple types of the vehicles, differences in financial accounting system of various corporations, and so on. Therefore the data has to be cleaned.Here we chose a certain range from 40 RMB to 60 RMB per load as the revenue of loaded movement between various freight terminals. The loading capacity of various type vehicles is allowed the certain range from 20 units to 30 units per vehicle. The cost of empty reposition is allowed the certain range from 200 RMB to 300 RMB per vehicle. We use the cost range from 2000 RMB to 3000 RMB per vehicles as purchase or lease cost of various vehicles.To evaluate the performance of the heuristics proposed, we generate randomly three sets of problem instances with different size of freight terminal and time periods, considering different number of vehicles types. Three sets of experiments are oriented to compare the performances of the hybrid simulated annealing algorithm with the CPLEX solver and VA-GHP approach.In the first set of contrast experiment, there are the heterogeneous fleet with three types of vehicles and five time periods. The number of the freight terminal is from 5 to 10. In the second set of contrast experiment, there are the heterogeneous fleet with five types of vehicles and ten time periods. The number of the freight terminal is from 11 to 20. In the third set of contrast experiment, there are the heterogeneous fleet with ten types of vehicles and thirty time periods. The number of the freight terminal is from 21 to 50.The test work has been done using a computer with Intel i5-5200CPU (2.2 GHz). The hybrid simulated annealing algorithm combining VA-GHP&USP-IVA has been coded in Microsoft Visual C++.


## 8.2. Performance Evaluation

In this section the computational experiments are performed using three methods to solve the MHFS-FSD model: HSAH program, GHP approach, and CPLEX solver. And we use two measures to evaluate the performance. The first one is the OPT, which is the value of the objective function obtained by the CPLEX solver, GHP approach, and HSAH program. OPT as one of the major criteria in assessing the performance of the HSAH, GHP, and CPLEX represents the profit generated by revenues for dispatching loaded vehicles, repositioning empty vehicles, and ownership costs for owning vehicle in planning horizon. The second measure of performance is the CUP time to run CPLEX solver, the GHP approach, and HSAH program for MHFS-FSD model. We aim at quantifying the benefits of the proposed method in two aspects: (1) the efficiency of the proposed method executed with no time limit imposed and (2) the efficiency of the proposed method executed with time limit imposed.(  1) Measure of Performance with No CPU Time Limit Imposed. In this section, we focus on the numerical studies on the performance of GHP approach, HSAH program and CPLEX solver executed with no CPU time limit imposed. We compare the performance of these three methods in two aspects: solution quality and computational efficiency. In Tables 1–3 we present the results obtained by solving randomly generated MHFS-FSD problems. We compare the solutions obtained by solving the MHFS-FSD problems using GHP approach, HSAH program, and CPLEX solver.Table 1
Performance for CPLEX, GHP, and HSAH applied to small size problem with no time limit imposed.
Problem size (freight terminal × vehicle types × time period) CPLEX GHP HSAH OS CPU time OFV CPU time RAT GAP OFV CPU time RAT GAP (RMB) (second) (RMB) (second) (%) (second) (RMB) (second) (%) (second) 5 × 3 × 5 219586 214.5 165894 13.6 24.45 200.9 199473 33.6 9.16 180.9 6 × 3 × 5 225465 235.7 184759 15.7 18.05 220.0 209484 37.5 7.09 198.2 7 × 3 × 5 230187 356.2 163028 13.9 29.18 342.3 190583 36.6 17.21 319.6 8 × 3 × 5 246586 248.3 190324 21.6 22.82 226.7 224756 45.6 8.85 202.7 9 × 3 × 5 250829 321.2 186206 25.8 25.76 295.4 208573 50.4 16.85 270.8 10 × 3 × 5 276585 314.5 201478 23.1 27.16 291.4 235859 55.6 14.72 258.9Table 2
Performance for CPLEX, GHP, and HSAH applied to medium size problem with no time limit imposed.
Problem size (freight terminal × vehicle types × time period) CPLEX GHP HSAH OS CPU time OFV CPU time RAT GAP OFV CPU time RAT GAP (RMB) (second) (RMB) (second) (%) (second) (RMB) (second) (%) (second) 11 × 5 × 10 727865 389.5 512034 89.7 29.65 299.8 708543 196.7 2.65 192.8 12 × 5 × 10 767289 397.2 536401 91.4 30.09 305.8 741478 183.4 3.36 213.8 13 × 5 × 10 876537 405.7 641963 103.4 26.76 302.3 835317 209.9 4.70 195.8 14 × 5 × 10 782685 412.3 541069 113.5 30.87 298.8 729294 220.3 6.82 192.0 15 × 5 × 10 813858 434.2 571648 73.8 29.76 360.4 787625 177.8 3.22 256.4 16 × 5 × 10 882385 451.8 642917 112.5 27.13 339.3 837674 229.6 5.07 222.2 17 × 5 × 10 926048 446.9 702541 100.9 24.14 346.0 892573 214.2 3.61 232.7 18 × 5 × 10 949874 467.1 731529 85.2 22.99 381.9 932074 192.1 1.87 275.0 19 × 5 × 10 998375 482.6 768359 78.4 23.04 404.2 948927 187.4 4.95 295.2 20 × 5 × 10 1010948 513.7 763542 104.3 24.47 409.4 978375 210.6 3.22 303.1Table 3
Performance for CPLEX, GHP, and HSAH applied to large size problem with no time limit imposed.
Problem size (freight terminal × vehicle types × time period) CPLEX GHP HSAH OS CPU time OFV CPU time RAT GAP OFV CPU time RAT GAP (RMB) (second) (RMB) (second) (%) (second) (RMB) (second) (%) (second) 21 × 10 × 30 3641457 1046.6 3162419 351.7 13.16 694.9 3592236 631.5 1.35 415.1 22 × 10 × 30 3960247 1153.7 3375864 363.2 14.76 790.5 3890437 654.2 1.76 499.5 23 × 10 × 30 4093746 1241.3 3392618 389.1 17.13 852.2 3962476 691.0 3.21 550.3 24 × 10 × 30 4172590 1306.9 3521607 412.5 15.60 894.4 4063582 719.6 2.61 587.3 25 × 10 × 30 4395799 1388.1 3862453 434.7 12.13 953.4 4298741 767.9 2.21 620.2 26 × 10 × 30 4673825 1449.7 3971641 443.8 15.02 1005.9 4495476 779.0 3.82 670.7 27 × 10 × 30 4691793 1533.8 4071358 503.4 13.22 1030.4 4586179 839.1 2.25 694.7 28 × 10 × 30 4762133 1592.4 4193062 536.7 11.95 1055.7 4703167 864.2 1.24 728.2 29 × 10 × 30 4792619 1640.3 4083017 576.3 14.81 1064.0 4698731 899.6 1.96 740.7 30 × 10 × 30 4821617 1731.6 4262843 598.6 11.59 1133.0 4803164 904.3 0.38 827.3 31 × 10 × 30 4865329 1796.1 4228019 602.6 13.10 1193.5 4829316 932.5 0.74 863.6 32 × 10 × 30 4887691 1874.2 4353762 647.9 10.92 1226.3 4796317 966.7 1.87 907.5 33 × 10 × 30 4963720 1924.0 4464907 712.4 10.05 1211.6 4850691 1043.2 2.28 880.8 34 × 10 × 30 5120379 1989.6 4526143 726.8 11.61 1262.8 5006479 1132.5 2.22 857.1 35 × 10 × 30 5436826 2046.5 4593011 779.3 15.52 1267.2 5267984 1164.6 3.11 881.9 36 × 10 × 30 5761844 2165.3 4612938 852.4 19.94 1312.9 5706914 1254.2 0.95 911.1 37 × 10 × 30 5839764 2221.2 4685714 837.0 19.76 1384.2 5784639 1275.6 0.94 945.6 38 × 10 × 30 6032875 2302.2 4703519 954.3 22.04 1347.9 5859217 1324.8 2.88 977.4 39 × 10 × 30 6361448 2399.5 4794808 989.4 24.63 1410.1 6263941 1396.1 1.53 1003.4 40 × 10 × 30 6673280 2556.4 4862714 1032.6 27.13 1523.8 6612476 1443.2 0.91 1113.2 41 × 10 × 30 6805218 2826.5 5182913 1058.1 23.84 1768.4 6569031 1504.1 3.47 1322.4 42 × 10 × 30 6961760 2973.3 5403917 1129.8 22.38 1843.5 6795213 1575.4 2.39 1397.9 43 × 10 × 30 7134752 3205.6 5632051 1186.7 21.06 2018.9 6948217 1614.5 2.61 1591.1 44 × 10 × 30 7446389 3360.5 5789830 1056.3 22.25 2304.2 7164206 1653.2 3.79 1707.3 45 × 10 × 30 7623796 3586.6 5913517 1219.4 22.43 2367.2 7501286 1741.8 1.61 1844.8 46 × 10 × 30 7853219 3888.2 6242683 1234.7 20.51 2653.5 7632689 1813.6 2.81 2074.6 47 × 10 × 30 7905256 4169.6 6590263 1196.4 16.63 2973.2 7693618 1856.8 2.68 2312.8 48 × 10 × 30 8134962 4275.3 6832657 1324.7 16.01 2950.6 8047603 1943.6 1.07 2331.7 49 × 10 × 30 8416732 4453.6 7207149 1717.2 14.37 2736.4 8267798 2066.2 1.77 2387.4 50 × 10 × 30 8601795 4679.1 7421678 1536.8 13.72 3142.3 8396419 2164.7 2.39 2514.4CPLEX solver is executed with no time limit imposed. The parameters used for the HSAH are as follows. The initial temperature is set byt0=2×F(ψ∗), where F(ψ∗) is the objective function value obtained by taking the initial solution solved with VA-IGP into objective function of MHFS-FSD model. The termination temperature is set to 1. The attenuation coefficient of temperature control parameter is assumed to be 0.9. The initial solution transformation frequency is allowed 80. The increasing rate of the solution transformation frequency is assumed to be 1.1.In order to facilitate problem analysis, we introduce two analytical indicators. The first denoted as RAT represents the deviation rate from the optimal solution. Let OFV be the objective function value with GHP approach and HSAH program. And let OS be optimal solution obtained with CPLEX solver. The RAT is calculated by(33)RAT=OS-OFVOSThe second is the CPU time gap denoted as GAP. It is calculated by the difference between CUP time of running the HSAH program or GHP approach and that of CPLEX solver. The CPU times of running HSAH program, GHP approach, and CPLEX solver are, respectively, denoted as CPUT.HSAH, CPU.GHP, and CPUT.CPLEX. The GAP is calculated by(34)GAP=CPUT.CPLEX-CPUT.HSAHGAP=CPUT.CPLEX-CPUT.GHPThe performance of small size problem is shown in Table1. Here we have 6 test problems. Each problem is, respectively, solved using the proposed GHP approach, HSAH program, and CPLEX solver. For the small size problem, the performance of GHP and HSAH is not very prominent when compared to CPLEX. And the performance of HSAH is better than GHP. The deviations of GHP from the optimal solution are the range from 18% to 30% while the deviations of HSAH from the optimal solution are the range from 7% to 20%. The CPU times to obtain the solutions are the range from 33.6 s to 55.6 s for HSAH, the range from 13.6 s to 23.1 s for GHP, and the range from 214.5 s to 356.2 s for CPLEX. The superiority of GHP and HSAH in computing time can be shown distinctly. The results of comparison are displayed in Figure 8.Figure 8
Results of comparison for small size problem.
(a)
Objective function value (b)
Deviation rate (c)
CPU running time (d)
CPU time gapThe performance of medium size problem is shown in Table2. Here we have 10 test problems. Each problem is respectively solved using the proposed GHP, HSAH, and CPLEX. For the medium size problem, the GHP and HSAH perform well when compared to CPLEX. The performance of HSAH is prominent than GHP. The deviations from the optimal solution are reduced to a reasonable range from 2% to 7% for HSAH and the deviations from the optimal solution range from 22% to 31% for GHP. One obvious observation is that the increase of CUP time to obtain the solutions using CPLEX begins to appear a trend of acceleration. But the rising step of CPU times to obtain the solutions using GHP and HSAH increases is relatively flat. So the superiority of GHP and HSAH in computation time is very obvious. The performance comparisons are shown in Figure 9.Figure 9
Results of comparison for medium size problem.
(a)
Objective function value (b)
Deviation rate (c)
CPU running time (d)
CPU time gapThe performance of large size problem is shown in Table3. Here we have 30 test problems. Each problem is respectively solved using the proposed GHP, HSAH, and CPLEX. For the large size problem, performance of solving the problem using the GHP and HSAH is more excellent than that of CPLEX. The deviations from the optimal solution for HSAH are limited to a very small range within 4% while the deviations from the optimal solution for the GHP are limited to range within 30%. But the CPU times for obtaining the solutions using CPLEX solver increase in a sharp speed due to the high memory requirements. The superiority of GHP and HSAH in computation time is particularly obvious. The performance comparisons are shown in Figure 10.Figure 10
Results of comparison for large size problem.
(a)
Objective function value (b)
Deviation rate (c)
CPU running time (d)
CPU time gap(  2) Measure of Performance with CPU Time Limit Imposed. In the experiment, we test the performance of the HSAH program and CPLEX solver with CPU time limit imposed. Since the performance of GHP approach is obviously inferior to that of HSAS program, GHP approach is no longer tested in the experiment. Here the HSAS program is fully executed. And let CPLEX solver be executed for the same CPU time. Then we compare the objective function value obtained by solving the MHFS-FSD problems using HSAH program and CPLEX solver. The tests are performed mainly in Tables 4–6.Table 4
Results for HSAH and CPLEX applied to small size problem with time limit imposed.
Problem size (freight terminal × vehicle types × time period) OFV (RMB) CPU time (second) CPLEX HSAH 5 × 3 × 5 53817 199473 33.6 6 × 3 × 5 57429 209484 37.5 7 × 3 × 5 61437 190583 36.6 8 × 3 × 5 68413 224756 45.6 9 × 3 × 5 73076 208573 50.4 10 × 3 × 5 78462 235859 55.6Table 5
Results for HSAH and CPLEX applied to medium size problem with time limit imposed.
Problem size (freight terminal × vehicle types × time period) OFV (RMB) CPU time (second) CPLEX HSAH 11 × 5 × 10 493276 708543 196.7 12 × 5 × 10 497625 741478 183.4 13 × 5 × 10 513410 835317 209.9 14 × 5 × 10 526683 729294 220.3 15 × 5 × 10 525914 787625 177.8 16 × 5 × 10 543728 837674 229.6 17 × 5 × 10 562807 892573 214.2 18 × 5 × 10 567436 932074 192.1 19 × 5 × 10 573614 948927 187.4 20 × 5 × 10 592460 978375 210.6Table 6
Results for HSAH and CPLEX applied to large size problem with time limit imposed.
Problem size (freight terminal × vehicle types × time period) OFV (RMB) CPU time (second) CPLEX HSAH 21 × 10 × 30 2561437 3592236 631.5 22 × 10 × 30 2746308 3890437 654.2 23 × 10 × 30 3012561 3962476 691.0 24 × 10 × 30 3348021 4063582 719.6 25 × 10 × 30 3585719 4298741 767.9 26 × 10 × 30 3602476 4495476 779.0 27 × 10 × 30 3680418 4586179 839.1 28 × 10 × 30 3764292 4703167 864.2 29 × 10 × 30 3791068 4698731 899.6 30 × 10 × 30 3822516 4803164 904.3 31 × 10 × 30 3849317 4829316 932.5 32 × 10 × 30 3886203 4796317 966.7 33 × 10 × 30 3965042 4850691 1043.2 34 × 10 × 30 4002837 5006479 1132.5 35 × 10 × 30 4125340 5267984 1164.6 36 × 10 × 30 4196205 5706914 1254.2 37 × 10 × 30 4273568 5784639 1275.6 38 × 10 × 30 4305217 5859217 1324.8 39 × 10 × 30 4362208 6263941 1396.1 40 × 10 × 30 4385426 6612476 1443.2 41 × 10 × 30 4412716 6569031 1504.1 42 × 10 × 30 4472930 6795213 1575.4 43 × 10 × 30 4502814 6948217 1614.5 44 × 10 × 30 4575356 7164206 1653.2 45 × 10 × 30 4598325 7501286 1741.8 46 × 10 × 30 4620819 7632689 1813.6 47 × 10 × 30 4673369 7693618 1856.8 48 × 10 × 30 4702923 8047603 1943.6 49 × 10 × 30 4785527 8267798 2066.2 50 × 10 × 30 4832509 8396419 2164.7From above three sets of contrast tests, the objective function values obtained by HSAH program are significantly better than those obtained by CPLEX solver with CPU time limit imposed.In sum, we can conclude that the HSAH is effective to obtain good solutions in relatively low run times for the MHFS-FSD. The results obtained indicate the benefits of using HSAH, helping MHFS-FSD to find better solutions in reasonable short CPU times, which are acceptable for practical applications.


## 9. Conclusions

The study described in this paper attempts to integrate vehicles allocation and fleet sizing decisions considering the vehicle routing of multiple vehicle types. We name the problem as multistage heterogeneous fleet scheduling with fleet sizing decisions (MHFS-FSD). In this paper, a mixed integer programming model has been developed for MHFS-FSD. As an NP problem, using traditional method for solving MHFS-FSD model is difficult and inefficient. A novel heuristic is considered to solve this problem.According to the characteristic of this problem, the matrix formulation denoting vehicle allocation scheme is provided. The approach of generating vehicle allocation scheme with greedy heuristic procedure (VA-GHP) as initial solution of problem is presented. On the basis of the initial solution generated by VA-GHP approach, the USP-IVA method to update the initial solution is developed. And then, incorporating VA-GHP and USP-IVA into simulated annealing algorithm, a novel heuristic called HSAH-GHP&IVA is proposed. Finally, we evaluate the performance of the hybrid simulated annealing heuristic in terms of traditional measure such as objective function and execution time. The results show that the proposed heuristic is able to obtain reasonably good solutions in short CPU times.Future research can focus on new heuristics to solve the MHFS-FSD. Another prospective research will be to extend the model of MHFS-FSD by introducing further realistic aspects and relaxing hard constraints.



---

*Source: 5364201-2019-01-10.xml*