Abstract
In this paper, a two-phase hybrid algorithm to address the problem of scheduling visits to customers and the vehicle routing problem of medication delivery to highly dependent patients is proposed. In the first phase, the issue of daily scheduling for a cluster of customers is solved via a flexible mathematical optimization model applied to different scenarios. The solution obtained in the first phase generates clusters of patients and a frequency routing problem that considers delivery periodicity, location, demand, service times, travel times, and daily load-balancing constraints. In the second phase, a hybrid metaheuristic approach including the synergy and constant iteration between simulated annealing and a record-to-record algorithm is applied to improve the initial solution obtained in the first phase. The effectiveness of the proposed algorithm is validated with real data obtained from a pilot project in Chile. The results are promising and demonstrate the efficiency of the proposed methodology.
1. Introduction
The problem considered in this paper is related to an important social issue: diseases that affect the independence of thousands of people when performing daily activities. Different studies conducted in several provinces in Chile have found that complex procedures at the time of medicament delivery cause high congestion levels in the hospital system and subject patients to diverse eventualities associated with their health conditions. New family compositions, high life expectancies, and social individualism place these patients in environments that are not prepared to meet their needs [1]. This study seeks to provide a solution to one of the most significant problems of this segment—home delivery of medications.
Progressive aging and an increase in the life expectancy of the elderly population, within the framework of new family compositions, require the adaptation of this population to challenging environments. A similar problem is experienced by people whose health conditions force them into permanent care situations. In Chile, 20% of the population has some degree of disability [1], and 8.3% of the population has severe impairment. These statistics indicate that approximately 1,082,965 people encounter performance issues in their daily lives. Given this context, only 35.6% of the population with disabilities receives assistance with daily activities, among which pick-up of medications from public healthcare centers is a primary problem [1].
The medication delivery problem is significant for all participants. This problem involves excessive waiting times that exceed four hours and the physical complexities of patients during the pick-up of medications, where the needs of patients are not always met. Partial solutions to the problem, such as medical home delivery services based on Internet requests, have been identified [2]. However, the lack of a logistics system, limited patient access to the Internet, and the combinatorial difficulty of the problem in Chile enable less than 100 monthly deliveries to a population with thousands of chronic and dependent patients.
In this paper, a heuristic algorithm is proposed for the problem of scheduling visits to customers and vehicle routing for the home delivery of medications and supply services. This algorithm aims to improve the quality of life of thousands of people with health conditions that make leaving their homes a significant challenge. The assignment of medications and supplies to drivers for subsequent delivery to patients is considered. The goal of the problem is to minimize delivery times and meet deadlines that are vital due to the criticality of the products (medications) to be delivered. The efficiency of the proposed methodology is tested with real data obtained from a pilot program in Chile.
The main contribution of this proposal is a management system that plans, coordinates, and controls each of the activities needed to perform a home delivery service for medicines and supplies for patients who are highly dependent, which is a highly complex problem. Patients may receive medication at frequencies that range from weekly and biweekly to monthly. Due to the severe dependence of the patients on these medications, the travel time between the locations of two patients and estimates of the average service time associated with each user are required, as medication delivery can require additional time that surpasses the travel time.
In general, the demand for medications and supplies per patient regarding weight and volume considering the size of an average truck (500 kg) is small; thus, the primary constraint is time, given that only 44 working hours per week are considered. The characteristics of the considered problem fit into a variant of the well-known vehicle routing problem (VRP) called the periodic vehicle routing problem (PVRP), which considers capacity and service time. Different from the well-known PVRP, the considered problem allows for exceeding the capacity of the vehicles with the payment of a penalty. Additionally, we consider the sum of the costs of overtime when the maximum duration of the routes is exceeded, the fixed costs of the days on which the visits are performed, and the costs related to the balance of days (days with similar workloads).
This paper is organized as follows. In Section 2, a literature review of similar problems related to scheduling visits to customers and vehicle routing problems for the delivery of medications is given. Section 3 describes the proposed methodology to address the considered problem. Section 4 and Section 5 provide the computational results and the conclusions and recommendations, respectively.
2. State of the Art
The considered problem is closely related to the well-known VRP, which is regarded as a nondeterministic polynomial-time NP-hard problem [3–5]. The operation costs for vehicle routing problems can be reduced between 10% and 20% using operational research techniques to facilitate planning [6–10]. The mileage can be reduced, staff and patient schedules can be optimized, established delivery deadlines can be satisfied, and the use of resources can be minimized.
The VRP is one of the most investigated combinatorial optimization problems [5, 11] due to its mathematical complexity, its importance in daily situations, and its relevance to any organization that provides home delivery of a product or service. The VRP can be considered a general category of a set of problems with a set of geographically dispersed customers for which the demand of products or services must be satisfied by a fleet of vehicles departing from one or more depots [12].
The solution of the VRP problem must specify the customers to be assigned to each vehicle and the order of the customers to be visited to minimize the total cost subject to a variety of constraints, such as vehicle capacity and dispatch times [5, 13, 14]. Similar related routing problems have been studied by Michalos et al., [15], Escobar et al., [9], Escobar et al., [10], and Bernal-Moyano et al. [16].
Different variants of the VRP incorporate various aspects relevant for industries or applications. For the considered problem, the formulation must exhibit specific characteristics, such as a planning horizon longer than one day, including service times, and a balanced generation of solutions. A PVRP, including some unique characteristics, can approximate the considered problem. The standard form of the PVRP has been extensively investigated by Cordeau et al. [17], Chu et al. [18], and Lacomme et al. [19] for a small number of customers. A substantial amount of data exists, which renders the solution of the problem by exact mathematical models very complex.
The PVRP is a generalization of the well-known VRP that attempts to determine an optimum set of daily routes for a given period. The PVRP provides the solution to two fundamental problems [20]: an assignment problem whose objective is the determination of a set of visiting days for each customer within a time horizon and a vehicle routing problem for each day to determine the best daily route. Each customer has a visit frequency during the time horizon, which implies that the number of solutions to the problem is large and exponentially increases with the number of nodes of the graph . The mathematical formulation of a PVRP must specify the level of service or visit frequency for each customer. Additionally, vehicle capacities and customer demands are required.
The majority of scientific articles published regarding the PVRP consider heuristic solutions. These types of approaches have been proposed by Beltrami and Bodin [21] and Russell and Igo [20]. Sophisticated algorithms to solve the PVRP have been introduced by Christofides and Beasley [22], Tan and Beasley [23], Gaudioso and Paletta [24], and Cordeau et al. [17]. Variants of the PVRP which consider time windows and service times (PVRPTW) have been extensively investigated by Pirkwieser and Raidl [25], Pirkwieser and Raidl [26], Pirkwieser and Raidl, [27], Michallet et al. [28], and Nguyen et al. [29].
From the perspective of the considered problem, the level of detail of the described formulation represents most of the required characteristics. However, in most studies, only problems with low periodicity are solved without considering the level of service and using a small number of nodes. Although no other research has yet addressed the problem described in this paper, some studies have addressed the problem of the home healthcare of patients, which is known as the home healthcare routing problem (HHCRP). Home healthcare (HHC) companies are widespread in European countries and some South American countries, such as Colombia. The main objective of the HHC problem is the home care of patients who are recovering from a disease or injury. Given the high transportation costs in the healthcare industry, research related to the optimization of logistics routes in home healthcare is essential. A state-of-the-art review of different problems associated with HHC is presented by Emiliano et al. [30] and Gutiérrez and Vidal [31].
A state-of-the-art review of the algorithms proposed for the HHCRP has been presented by Fikar and Hirsch [32] and Cissé et al. [33]. Fikar and Hirsch [32] present a general approach to the different routing and programming methodologies based on the problem settings. Modern optimization techniques applied to the problem are highlighted, and some future studies are discussed. Cissé et al. [33] present different operational research models for the home healthcare routing and scheduling problem (HHCRSP). The focus of this study is to show the formulation of various objective functions and constraints. The HHCRP has been addressed using different exact and approximate methodologies. Exact deterministic algorithms for the solution of the HHCRP have been proposed by En-nahli et al. [34], Decerle et al. [35], and Manerba and Mansini [36]. Exact stochastic algorithms have been proposed by Shi et al. [37]. The proposed exact approaches for the HHCRP have been able to solve only small-size instances (up to 50 customers).
Trajectory-based metaheuristics have been proposed by Steeg and Schröder [38], Trautsamwieser and Hirsch [39], Liu et al. [40], and En-nahli et al. [41]. Steeg and Schröder [38] propose a solution scheme based on a combination of the strengths of constraint programming and a variable neighborhood search (VNS) metaheuristic method. The problem considers the minimization of the number of nurses that visit a patient during the planning horizon. Trautsamwieser and Hirsch [39] suggest a mathematical model for the HHCRP and a metaheuristic based on a variable neighborhood search. The objective function considered minimizes the travel times of nurses. A tabu search considering different feasible and infeasible neighborhoods has been proposed by Liu et al. [40]. In this study, the results obtained from an HHC company reveal that the proposed scheme reduces the total cost and optimizes the vehicle workload balance. En-nahli et al. [41] propose a metaheuristic method based on an iterated local search combined with a random variable neighborhood search. Some disadvantages of these approaches include slow convergences rates, being trapped in local optima due to the quality of the initial solution and the use of a large number of parameters in the approaches, making the experimental design for the tests difficult.
Population algorithms for the HHCRP have been proposed by Akjiratikarl et al. [42], Shi et al. [43], and Decerle et al. [44]. Akjiratikarl et al. [42] present a novel application of the particle swarm optimization (PSO) metaheuristic for the scheduling problem of home care workers. The objective function considered minimizes the distance travelled while satisfying the capacity and time window constraints [14, 45]. A stochastic version of the problem using a hybrid genetic algorithm with simulation is proposed by Shi et al. [43]. This version of the problem considers the service time and the travel time as probabilistic parameters. Decerle et al. [44] introduce a memetic algorithm to evaluate a multiobjective approximation scheme using several instances and support decision-making. Some of these proposed approaches have drawbacks, such as the use of many complex operators [42], the use of many parameters to be tuned [44], premature convergence [44], difficulties for convergence [42], some unpredictable results [43], and long convergence times [43]..
Approaches combining exact and approximate techniques have been developed by Bertels and Fahle [46], Liu et al. [47], and Allaoua et al. [48].
3. Materials and Methods
The considered problem, which is a generalization of the PVRP, can be solved by using mathematical models or heuristic approaches. Exact methods can optimally solve only small- or medium-sized instances. On the other hand, heuristic algorithms can solve large problems, such as the considered problem.
The decisions involving our problem are as follows: (a) the assignment of patients to visit days, (b) the assignment of patients being visited to a particular vehicle each day, and (c) the sequencing of patients assigned to a particular day and vehicle. Solving the three problems together is very complicated. The problem is solved by using a two-phase algorithm, in which each phase independently considers a subproblem. In the first phase (long-term decisions), problem “(a)” is solved considering the days to visit the patients. Second phase (short-term decisions) solves problems “(b)” + “(c)” by considering decisions related to the assignment of the vehicles and the order of visiting the patients.
For Phase 1, the solution of the problem of scheduling visits to customers in clusters for each day of the horizon, according to the location, delivery frequency, and service time, is achieved by separating the original problem into smaller subproblems. This phase is solved exactly using an integer linear programming problem.
Regarding Phase 2, the VRP with capacity and service time is solved for each of the daily clusters generated in Phase 1. This phase uses a heuristic based on the techniques proposed by Groër et al. [49].
Phase 1 (solution of the visit scheduling problem). The process of clustering the patients is performed according to the delivery frequency, location, service time, and delivery time using an integer linear programming mathematical model.
Sets. See Table 1.
Constants. See Table 2.
Parameters. See Table 3.
The visits to the patients are always performed on the same days of the week (fixed days). For example, a customer with weekly visits who is visited on the first Monday of the month must be visited on the Monday of every week. If a customer has a fortnightly frequency and is visited on the Wednesday of the second week, this customer must be visited on the Wednesdays of the second and fourth weeks of each month.
The number of visits to perform in a month is defined by the function . The function normalizes the current day to the day that the current represents (first day of visit) for the customer (the customer is always visited on the same day of the week). For example, if a customer has an assignment of weekly frequency, day 1 (Monday of the first week) represents days 6 (Monday of the second week), 11 (Monday of the third week), and 16 (Monday of the fourth week). If a customer has an assignment with biweekly visits, day 4 (Thursday of the first week) represents day 14 (Thursday of the third week).
Decision Binary Variables
Decision Variables : Lower and upper bounds in radians for day . : Lower and upper bounds for daily deliveries.
Objective Function
The objective function includes the following information:(a)Cost of extra weight for each day (seeks to avoid exceeding the load capacity of each vehicle).(b)Cost of the extra hours required in the medication delivery process. The cost is considered to be the maximum number of hours allowed by law in one day. A penalty is incurred if delivery is performed outside regular service hours.(c)Cost of worked days. A driver’s agenda can be compacted using the maximum hours available in a day and generating days off. The cost is equivalent to a typical working day.(d)Cost of balancing orders. This function enables approximately the same number of trips to be performed daily, is conditional, and employs large clusters of patients, similar to the problem under consideration. In particular, the proposed model considers the difference between the days with the largest and the lowest number of deliveries . This result corresponds to the “penalization term,” which is multiplied by the cost of balancing the daily work (). The cost of balancing orders is based on the importance given to unbalanced routes concerning other costs.(e)Cost of clustering. The general idea is to cluster the customers who are near each other using an angular distance clustering constraint.
Constraints(a) Belonging to each frequency: each customer on day must have an associated delivery frequency: (b) Frequency clustering: each customer who has a delivery on day will have the next delivery on the same day of week , which corresponds to the delivery frequency. The assignment will be performed by the frequency normalization function :(c) Vehicle capacity constraint: the daily capacity () of each vehicle should not be exceeded; the daily capacity is defined as the sum of the demands () for customers on day . When extra capacity is required, a fee is established for buying additional capacity ().(d) Time constraint: working hours () should not be exceeded; the working hours are defined as the sum of the service time () for customer and the mean time estimator (). When additional time is required, a daily fee for paying extra hours is established (). (e) Daily clustering constraint: when the problem is small, customers can be clustered into the smallest number of days; this constraint generates clusters of customers and maximizes the use of the daily working hours. The number of working days is minimized. is defined as the cardinality of the set of customers . (f) Balancing constraint: when the problem is large, the workload of the drivers must be balanced by dividing the number of deliveries that must be performed per day into equivalent quotas. The lower limit is dynamically defined by , and the upper limit is dynamically defined by . (g) Distance clustering constraint: the centroid is established as the average of all positions () for each customer . Then, the angle between the position difference of customer and the centroid is calculated in radians. The angles are limited to the left by a lower bound () and limited to the right by an upper bound (). (h) Variables domain: the variables are binary, continuous, and greater than or equal to zero.
Phase 2 (solution of the vehicle routing problem). In this phase, certain amounts of vehicle routing subproblems are obtained due to the periodicity of the general problem. Each customer has demands for products and service times.
A procedure based on the VRPH software library [49] containing local and global search heuristics and metaheuristics for vehicle routing problems is employed to solve the problem. This method uses a variant of the Clarke and Wright [50] heuristic for determining an initial solution (Escobar et al., 2003), applying the simulated annealing metaheuristic [51] and the record-to-record metaheuristic [52] and iteratively obtaining a solution within an acceptable computational time.
The proposed algorithm combines the metaheuristics described in a heuristic procedure adapted from Escobar et al. [53]. The programmed routines are as follows: (i)vrp_initial: this routine uses a variant of the Clarke-Wright algorithm(ii)vrp_record_to_record: this routine is an implementation of record-to-record travel(iii)vrp_simulated_annealing: this routine is an implementation of simulated annealing A feasible basic solution is created with the vrp_initial routine, and the vrp_simulated_annealing and vrp_record_to_record routines are iteratively applied until no improvements are achieved. The procedure is executed in several parts as a two-phase hybrid algorithm that generates significant improvements for a depot. The procedure is described in Algorithm 1.
| ||||||||||||||||||||||||||
The main differences between the approach used and the algorithm presented in Escobar et al. [53] for the CLRP are as follows: (i) the use of the complete graph instead of the granular search space, (ii) the consideration of the service times within the procedure, (iii) the possibility of considering only feasible solutions, and (iv) the use of the proposed procedure as the main solution for the heuristic part (Phase 2) of the former algorithm.
Figure 1 shows a general diagram of the proposed methodology, where Phase 1 is composed of mixed integer linear programming exactly solved by Cplex. The output of Phase 1 is the variables of the model, which indicate whether customer is visited on day . Phase 2 executes a heuristic algorithm with the assigned customers for each day. This approach generates an initial solution via the Clarke-Wright algorithm, which is later improved via a simulated annealing and a record-to-record approach until no improvements are found. Finally, the routes are obtained for each day.
4. Computational Results and Discussion
The proposed algorithm is programmed in C++. For the solution of Phase 1 (mathematical model), an IBM ILOG Concert Technology library with Cplex 12.6.2.0 is employed. The implementation of Phase 2 is performed entirely in C++ with a release adapted from the library proposed by Escobar et al. [53]. From a hardware perspective, the instances are solved using a MacBook Air Mid 2012 with a 1.8 GHz Intel i5-3427U processor and 4 GB of RAM.
4.1. Implementation Details
The instances are generated using a pilot program in Chile for the real problem, which contains location, demand, service time and medication, and supply delivery frequency data for 800 patients. Two types of instances are generated to obtain statistically representative data from the database. In the first instance, for the problem of scheduling visits to patients (Phase 1), modifications to the design of Cordeau et al. [17] are performed. In Phase 2, the configuration proposed by Christofides and Beasley [22] is applied to the vehicle routing problem.
The service time () ranges between three minutes and twelve minutes, depending on the complications that a patient may encounter, such as opening doors, finding documents, and receiving instructions. For the estimation of service times, a regression analysis is performed using the data obtained from the real case. However, the service time is probabilistic. Thus, a regression function based on the medication demand (), delivery frequency (), and patient age () is initially employed, since these known variables can have a strong impact on the service times. The demand depends on the number of medications and supplies that have to be delivered and validated; the higher the frequency is, the higher the probability that the patient has a severely disabling disease, and the patient age is correlated with increasing difficulties for a patient to move or find documentation. From a total of 800 patients, 58 observations are investigated, and a multiple correlation coefficient of 0.872 and an determination coefficient of 0.761 are obtained, which generate the coefficients listed in Table 4.
The time function is estimated as follows: For the calculation of distances, the Manhattan distance or taxicab distance [54] is employed. In particular, the Manhattan distance has benefits in well-designed cities because this method approximates the real distance better than the Euclidean distance. For the travel time calculations, a constant average speed of 35 [km/h] is assumed.
For the transformation from latitude and longitude in the geographic system to the Cartesian coordinate system, the Earth is assumed to be completely spherical, and the depot is established as the origin (0,0). Thus, the distance measurement is relative to this point. Let and be a pair of geographic coordinates measured in radians, defined by latitude and longitude, and let R be the radius of Earth in meters. The transformation into Cartesian coordinates is specified as follows:In the mathematical model, some parameters that determine the cost or specific weight to perform a specific action are considered. These parameters have to be adequately defined to guide the optimization according to the priorities and requirements of the problem. These parameters are specified with a specific weight in increasing order according to their importance. Table 5 summarizes the parameters used in the mathematical model.
4.2. Experimental Design
Two initial experiments are performed with an instance adapted from the literature (based on the M-n101-k10 instance proposed by Christofides and Beasley, [22]) and a randomly generated instance (with clusters of customers located in areas that are far apart, which enables graphic visualization of the cluster assignment for each day) before the tests are performed with a set of real data. Both initial instances have 100 customers. The first instance has constant service times and delivery frequencies. The second instance considers random customers, demands, and frequencies. The other clustering and balancing constraints are deactivated. Two cases are evaluated for each experiment.
(a) Experiment 1: Clustered Instance. In this experiment, the customers have been previously clustered to test constraints (9). Two cases are evaluated: the first case does not have the proximity clustering constraint (Figure 2(a)), and the second case employs this constraint (Figure 2(b)). Figure 2 shows the spatial location of each customer , where each day of the week is represented as follows: L = Monday, M = Tuesday, X = Wednesday, J = Thursday, and V = Friday.
(a)
(b)
In this experiment, the value of the objective function for the case without proximity clustering is the maximum value due to complete disorder in the clustering of the points. For the case with clustering, these points are ordered in the best possible manner, as graphically depicted.
(b) Experiment 2: Random Instance. In this experiment, a completely random instance, in which the points are widely dispersed throughout a plane, is evaluated; the demand and periodicity of the orders are random. Figure 3 shows the spatial location of each customer , where each day of the week is represented as follows: L = Monday, M = Tuesday, X = Wednesday, J = Thursday, and V = Friday.
(a)
(b)
In this experiment, the objective function value for the case without proximity clustering (Figure 3(a)) is high due to the high entropy in the grouping of the points. In the case with proximity clustering (Figure 3(b)), these points are ordered and produce a significantly lower objective function value than the previous case with a high computational time.
In small- and medium-sized instances, the use of the daily clustering constraint is more advisable than the application of the time-balancing constraint, and the use of both constraints is recommended. Daily clustering is conditional, since the results are only favorable for some instances, in which the orders can be fulfilled in fewer days. In other cases, daily clustering unbalances the problem when the use of fewer working days is favored.
Daily clustering is convenient for small- and medium-sized instances, since all orders can be distributed within fewer days while complying with the periodicity constraints. The measurement of the objective function value is performed according to the cost of working one day. The computational time increases compared with the other constraints, although this increase is negligible.
4.3. Results Obtained for a Real Instance
The proposed model is executed with a real instance of 800 customers, among which 32 customers (4.00%) receive weekly care, 78 customers (9.75%) receive biweekly care, and 690 (86.25%) customers receive monthly care. A total of 974 deliveries and 1,014 trips must be performed. The cumulative service time is 144.99 hours, and the cumulative demand is 851.51 kg per month. The number of working days per month is 20, the number of working hours is nine, eight hours are considered for delivery and loading, and one hour is considered for rest. A maximum of 2.4 extra hours can be allocated if required. This finding implies a total of 160 hours for delivery and a maximum of 48 additional hours. The vehicle has a daily capacity of 650 kg and an unrestricted schedule.
In the first instance, a mathematical model with all the constraints except for daily clustering is employed, since this model is applied to smaller problems. Due to a large number of possible combinations and the low-quality lower bound, obtaining optimum solutions within an acceptable period is difficult. However, a high-quality, feasible solution is obtained in the computational time of three hours with a gap of 25.03%. The objective function value corresponds to 20,144.26, and the lower bound value is 15,090.85. Table 6 shows the size of the problem, that is, the numbers of variables and constraints.
The heuristic algorithm is initially validated in two instances using service times. The first instance has 41 customers, which is similar to the number of daily patients in the real problem; the second instance is adapted from the study of Groër [55], with 500 customers. The initial results obtained with the Clarke-Wright heuristic are compared with the results obtained with the proposed hybrid procedure. The solution for the 800-customer instance is pursued.
From Phase 1, 20 routing problems that correspond to each scheduling day are obtained; each problem considers the capacity, demand, and service time for each customer. Initially, a solution is obtained using the Clarke-Wright algorithm. Then, the algorithm is improved using a two-phase algorithm that alternates between simulated annealing and record-to-record metaheuristics until a near-optimal solution is obtained. The average execution time of the algorithm in the second phase per day is 0.43 seconds. Regarding the real case, the algorithm is capable of decreasing the travel time by 59.5% and reducing the total working time by 8.5%.
The proposed approach allows solving the problem within short computing time (about 4 hours) in comparison with a real planning time by using a conventional method which usually considers more than three days. Besides, the proposed approach contemplates the visit of all the patients by considering their attention day, which is not possible with a manual planning process. Finally, Figures 4 and 5 show the behavior of the performance of the first and the second phase as a function of the number of nodes (first phase) and as a function of the computing time (second phase), respectively.
5. Concluding Remarks
In this article, a two-phase heuristic algorithm is proposed for the problem of scheduling visits to customers and vehicle routing for the delivery of medications to patients with a high degree of dependence.
In the first phase, the problem of grouping customers into clusters is solved via an integer linear programming model that generates groups of patients for each day using criteria such as the periodicity, distance, and travel and service times. The model is satisfactory for instances of any size when frequency clustering, balancing, and daily clustering constraints are employed. In the second phase, a hybrid metaheuristic that includes the synergy and continuous iteration of algorithms based on simulated annealing and record-to-record travel is proposed. This algorithm obtains near-optimal solutions in times of less than three seconds per day.
Future research is proposed to integrate this methodology with inventory decisions, avoid out-of-stock scenarios, and reduce the periodicity of weekly deliveries of medications to a higher number of patients.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work has been partially supported by project CONICYT FONDECYT 11150370, “Grupo de Logística y Transporte,” at Universidad del Bío-Bío and Universidad del Valle, Colombia. This support is gratefully acknowledged.