Integrated Harvest and Distribution Scheduling with Time Windows of Perishable Agri-Products in One-Belt and One-Road Context

The unique characteristics of perishable agri-products are a short lifespan and rapid quality deterioration. This establishes the need to significantly reduce the time from harvest to distribution. These features require reducing the processing time from harvest to distribution to being as short as possible. In this study, we focus on an integrated perishable agri-products scheduling problem that combines harvest and distribution simultaneously, with the purpose of reducing processing time and quality decay. We propose this problem as a mixed integer nonlinear programming model (MINLP) to optimize the harvest time and the vehicle routing to consumers, and this MINIP is formulated as a vehicle routing problem with time windows (VRPTW). We introduce a big M method to transform the nonlinear model into a linear model, then apply CPLEX to solve the transformed model. Numerical experiments and sensitive analysis are conducted to verify the efficiency of the proposed model and to provide managerial insights.


Introduction
In recent years, the Chinese government has adopted a number of policies to promote the opening up of the country and economic cooperation with other countries under the One-Belt and One-Road (OBOR) national strategy.This strategy can not only accelerate sustainable development but can also enhance international trade.China is one of the largest agricultural producing countries in the world.According to reports, the export of agri-products was at ¥11,209.1 billion in 2016, which was an increment of up from ¥503.49 billion as compared to the exports in 2015 [1].Perishable agri-products are not only necessary for people's daily lives but also an essential source of increasing income for peasant farmers.However, the short lifespan and quality decay nature of perishable agri-products make its deterioration rapid over time.According to an FAO report, the decay rate of perishable agri-products is about 1.7-5% in developed countries, while it is up to 20-30% along the farm-to-door supply chain in China [2].Moreover, the ratio of social logistics costs in the gross domestic product is almost two-times higher than developed countries [3].Therefore, reducing quality decay and shortening the farm-to-door supply chain so as to maintain freshness and reduce the logistics costs has been an essential problem for perishable agri-products trading in the OBOR context.
The integrated perishable agri-products harvest and distribution scheduling provides a useful way to compress the farm-to-door supply chain.The decision problems shed light on harvest time and distribution routes in the farm-to-door supply chain according to the harvest location, the consumer geographical location and demand time windows, and the quality deterioration over time.In this study, we aim to answer two questions: (1) How to optimize the harvest time and its relationship with distribution time and the arriving time to consumers; (2) How to make integrated perishable agri-products harvest and distribution scheduling with time windows.
Our contributions in this study include three aspects: (1) present a quantitative method to model the quality decay process of perishable agri-products; (2) propose a mixed integer nonlinear programming model to formulate the integrated harvest and distribution scheduling with time windows; (3) by adding auxiliary decision variables, the proposed nonlinear model is transformed into a linear model and solved by CPLEX, and carries out numerical experiments to verify the validity and efficiency of the proposed model.
The remainder of this paper is organized as follows.Section 2 reviews relevant literature.Section 3 proposes a mixed integer nonlinear programming model.Section 4 reports the numerical results.Finally, Section 5 summarizes the concluding remarks and future works.

Literature Review
Our study is related to three streams of research literature, namely, the quantitative measure of agri-products quality decay, agri-products harvest optimization, and agri-products distribution optimization.
In the research stream on the quantitative measure of agri-products quality decay, the proposed methods can be grouped as two aspects, i.e., the exponential function and the linear function.For the exponential function, Giri et al. [4] considered exponential quality decay in the supply chain with multiple manufacturers to decide the optimal pricing strategy.Ghare [5] proposed a negative exponential function relationship between quality loss and time variation.Rong et al. [6] stated that the meat quality deterioration decreases exponentially with time.Wang and Li [7] addressed perishable food pricing strategies and proposed an exponential quality decay weighting function.On the other hand, the linear quality decay function has been adopted popularly in the literature.Mohapatra and Roy [8] considered a uniform quality loss of perishable products over time.De Keizer et al. [9] addressed an average quality decay function for heterogeneous products.Nekvapil et al. [10] introduced a Raman coefficient into a linear decreasing quality decay to assess the freshness of different citrus groups.Besides, some other quantitative measures for agri-products quality decay were also investigated, for example, the Tobit regression model in Chen et al. [11], Weibull distribution in Qin et al. [12], the weight analysis method in Kasso et al. [13], and the Gompertz model and Arrhenius equation in Bruckner et al. [14].
For the agri-product harvest optimization problem, Thuankaewsing et al. [15] investigated the harvest scheduling problem and proposed an ANN-based sugarcane prediction model to predict the sugarcane harvest time.Higgins et al. [16] pointed out that optimization on harvest date could increase sugars' productivity and profitability.Lamsal et al. [17] addressed the integrated sugarcane harvest and delivery scheduling problem.Caixeta-Filho [18] studied the optimal orange harvest time.Supsomboon and Niemsaku [19] considered an integrated cultivation and harvest planning problem.Ferrer et al. [20] and Arnaout and Maatouk [21] considered the grape quality loss function in grape harvest planning.Jena and Poggi [22] studied optimal harvest time for maturation peaks.Accorsi et al. [23] showed that weather condition control for cherry harvest can reduce energy cost in the cold chain distribution.For more details related to agri-product harvest optimization, please refer to the literature review by Kusumastuti et al. [24].
In the research stream of agri-products distribution optimization, some extended works focused on the vehicle routing problem with time windows (VRPTW) of agri-products.Mohammed et al. [25] investigated the vehicle routing problem (VRP) of perishable food.Wang et al. [26] addressed a VRPTW for perishable products, with the objectives of minimizing distribution cost and maximizing freshness.Zhang and Chen [27] proposed a solution to the VRPTW in multi-product frozen food distribution considering unit volume, price, and perishable coefficient.Hsu and Chen [28] focused on analyzing a joint food distribution system operation strategy to the VRPTW of food delivery.Some academic work has covered methods for producing an optimal model of agri-products distribution, including Song and Ko [29], who considered various types of vehicles to distribute perishable products and formulated this problem as a mixed integer nonlinear programming (MINLP) model.Rocco and Morabitol [30] presented a linear programming model to optimize the logistic planning for the Brazilian tomato processing industry.Ahumada et al. [31] proposed a two-stage stochastic programming model for fresh agri-products production and distribution scheduling.Amorim and Almada-Lobo [32] formulated a perishable food delivery problem as a multi-objective programming model, with the objective of minimizing distribution costs and maximizing freshness.Lamsal et al. [33] formulated a perishable agri-products distribution problem as a mixed integer programming model and proposed a two-phase algorithm to solve the model.
There are many studies addressing the integrated optimization problems in the area of agri-products operations.Ahumada and Villalobos [34] established an operational model of harvest and distribution of fresh agri-products and concluded that increasing vehicles could reduce the freshness cost.Rong et al. [13] presented a mixed integer linear model for food production and distribution planning in a food supply chain.Soto-Silva et al. [35] proposed a programming model to optimize the purchasing, transporting, and storing decisions for a fresh apple processing plant.Wu et al. [36] and Chen et al. [37] addressed the perishable food production and delivery problem, and Wu et al. [38] pointed out that the optimal arrangement of vehicles can reduce the loss of perishable food significantly.
In summary, incorporating quality decay into agri-products logistics system has attracted more and more attention in the area of agri-products operations management.The existing studies mainly focused on the agri-products logistics optimization in the distribution phase, and few studies have investigated the integrated scheduling problem that combines the farm side with the consumption side in agri-products logistics systems.However, with the increasing awareness of freshness and quality by consumers, how to reduce the perishable agri-products quality decay through breaking the barriers in the farm-to-door supply chain has become more and more important.Therefore, this study aims to endeavor to bridge a gap in literature by integrating harvest and distribution scheduling together and optimizing the harvest time and vehicle routing for perishable agri-products.

Model Formulation
We considered an integrated harvest and distribution scheduling with time windows for perishable agri-products and modeled this problem as a VRPTW.The objective was to minimize the total costs, including the average quality decay cost in the harvest phase, the average quality decay cost in the distribution phase, and the distribution cost.

Problem Statement
The logistics network in this study consisted of one farmer, multiple consumers located in different places, and multiple homogeneous vehicles.It was given that the farmer had only one harvest location and provided one kind of perishable agri-product.The harvest location had enough capacity to meet all consumers' demand.Moreover, we assumed each consumer was served only by a single vehicle.All vehicles left and returned to the center in the harvest location, and the number of used vehicles could not exceed that which the center had available.Besides, in the distribution phase, we only considered the harvest time and distribution time, and neglected the loading and unloading time.
Let t k i denote the time of vehicle k arriving at consumers i, τ ij denote the distribution time from node i to node j, t s k and t f k denote the starting harvest time and the ending harvest time of perishable agri-products distributed by vehicle k respectively, and t k 0 denote the starting distribution time of vehicle k leaving at the harvest location.In this study, we assumed the time gap between the ending harvest time and the starting distribution time, i.e., t f k was equal to t k 0 .Figure 1 demonstrates the harvest and distribution scheduling process.Vehicle k served consumers 1, 2, and 3 starting from harvest location 0. The optimal decision problems focused on what time to do harvest operation, what time the vehicle leaves the harvest location, how many vehicles needed to serve the consumers, and how many consumers are served by the same vehicle and at what optimal sequence.In this study, we aimed to capture and formulate these decision questions.
Note that, for the explanation of parameters and variables used in the proposed model, please refer to Appendix A. serve the consumers, and how many consumers are served by the same vehicle and at what optimal sequence.In this study, we aimed to capture and formulate these decision questions.Note that, for the explanation of parameters and variables used in the proposed model, please refer to Appendix A.

Objective Function
We considered three different objectives in this study, i.e., the average quality decay cost in the harvest phase, the average quality decay cost in the distribution phase, and the distribution cost.

Average Quality Decay Cost in the Harvest Phase
Figure 2 depicts the quantity variation from to , in which denotes the demand of consumer and denotes the shadowed area.In this study, we considered the quality decay of perishable agri-products from the starting harvest time.Let γ denote the quality decay rate.Inspired by the study in [35], we can calculate the average quality decay from the starting harvest time to the ending harvest time as shown in Figure 2. The average quality decay can be calculated as .To optimize the quality decay cost during the harvest time period, the key problem focused on how many consumers and their delivery sequence were served by a same vehicle.Therefore, in the harvest phase, the optimal decisions were how much total time was needed to harvest and when to start the harvest operation, so as to meet the time window and demand requirement for each served consumer in the following distribution phase.
Let denote the average quality decay cost in the harvest phase.The calculation is given in Equation (1).

Objective Function
We considered three different objectives in this study, i.e., the average quality decay cost in the harvest phase, the average quality decay cost in the distribution phase, and the distribution cost.

Average Quality Decay Cost in the Harvest Phase
Figure 2 depicts the quantity variation from t s k to t f k , in which d i denotes the demand of consumer i and s denotes the shadowed area.In this study, we considered the quality decay of perishable agri-products from the starting harvest time.Let γ denote the quality decay rate.Inspired by the study in [35], we can calculate the average quality decay from the starting harvest time t s k to the ending harvest time t f k as shown in Figure 2. The average quality decay can be calculated as sγ serve the consumers, and how many consumers are served by the same vehicle and at what optimal sequence.In this study, we aimed to capture and formulate these decision questions.Note that, for the explanation of parameters and variables used in the proposed model, please refer to Appendix A.

Objective Function
We considered three different objectives in this study, i.e., the average quality decay cost in the harvest phase, the average quality decay cost in the distribution phase, and the distribution cost.

Average Quality Decay Cost in the Harvest Phase
Figure 2 depicts the quantity variation from to , in which denotes the demand of consumer and denotes the shadowed area.In this study, we considered the quality decay of perishable agri-products from the starting harvest time.Let γ denote the quality decay rate.Inspired by the study in [35], we can calculate the average quality decay from the starting harvest time to the ending harvest time as shown in Figure 2. The average quality decay can be calculated as .To optimize the quality decay cost during the harvest time period, the key problem focused on how many consumers and their delivery sequence were served by a same vehicle.Therefore, in the harvest phase, the optimal decisions were how much total time was needed to harvest and when to start the harvest operation, so as to meet the time window and demand requirement for each served consumer in the following distribution phase.
Let denote the average quality decay cost in the harvest phase.The calculation is given in Equation (1).To optimize the quality decay cost during the harvest time period, the key problem focused on how many consumers and their delivery sequence were served by a same vehicle.Therefore, in the harvest phase, the optimal decisions were how much total time was needed to harvest and when to start the harvest operation, so as to meet the time window and demand requirement for each served consumer in the following distribution phase.
Let obj 1 denote the average quality decay cost in the harvest phase.The calculation is given in Equation (1).
where t f k and t s k are continuous decision variables and x k i is a binary decision variable, which equals 1 if the consumer i is served by vehicle k and 0 otherwise.Therefore, obj 1 is a nonlinear objective function which could hardly be solved.We adopted the big M method to transform it into a linear model by adding auxiliary decision variables.The big M method is an effective and regular approach used to deal with the nonlinear integer optimization and has been widely adopted to convert nonlinear mixed integer decision variables into a linear one [38].The formulated function is given in Equation (2). where They are subject to the feasible constraints given in Equation (3).

Average Quality Decay Cost in the Distribution Phase
In this objective, the quality decay cost was calculated from the vehicle leaving time at the harvest location 0 to the arriving time at consumers i.Let obj 2 denote the average quality decay cost in the distribution phase, which can be calculated in Equation (4).
Similarly, obj 2 is also a nonlinear objective function which can be transformed into an equivalent linear form by using the big M method.Let p k i = t k i x k i and l k 0i = t k 0 x k i .obj 2 can be rewritten in Equation ( 5), with the feasible constraints given in Equation ( 6).

Distribution Cost
This objective contained fixed vehicle cost and distribution cost.Here, λ k ij is a parameter, and if node i = 0, then λ k ij = λ + α, otherwise λ k ij = λ.Let λ denote the unit distribution cost of vehicle k, α denote the fixed cost of vehicle k. x k ij is a binary decision variable, taking a value 1 if the arc (i, j) belongs to the route of vehicle k and 0 otherwise.The distribution cost can be calculated in Equation (7).

Constraints
We formulated the model as an VRPTW problem and constraints are given as follows.
Constraints ( 8)-( 12) were constraints of VRP.Constraints ( 8) and ( 9) limited each consumer to be visited by only one vehicle.Constraint (10) was a network flow conservation constraint.Constraints (11) and ( 12) limited each vehicle departing and returning the harvest location to at most one time.Constraint (13) referred to the vehicle k that will serve consumer i if it is a node on the corresponding route.Constraint ( 14) was a constraint of vehicle capacity.Constraint (15) was a constraint of vehicle quantity.Constraint (16) defined the relationship between the starting harvest time and the ending harvest time for vehicle k.Constraints ( 17) and (18) denoted the relationship between finishing harvest time, starting distribution time, and arriving time at consumer i.Constraint (19) limited the time of vehicle k arriving at consumer i.Note that, constraint (19) can eliminate the case of subtours because, 20) and ( 21) were constraints on time windows at harvest location and each consumer.Constraint ( 22) limited x k ij , x k i as binary decision variables.

Programming Model
According to the discussions in Section 3.2 and 3.3, the total objective function was to minimize obj 1 , obj 2 , and obj 3 by subjecting them to constraints (3), ( 6), ( 8)-( 22).However, since the magnitude of obj 1 , obj 2 , and obj 3 were different, we adopted the Min-Max normalization method to normalize the three objectives in the same magnitude, i.e., obj * i = Therefore, the integrated harvest and distribution scheduling model after normalization can be written as the following: (3), ( 6), ( 8)-( 22).
Note that model ( 23) is a VRPTW problem, which has proved to be an NP-hard problem [39].If constraints (19), (20), and ( 21) are removed, model ( 23) is equivalent to the VRP problem.Moreover, if let e 0 = e i = 0, l 0 = l i = M, and only one vehicle used, model ( 23) is equivalent to the travelling salesman problem (TSP).If the vehicle capacity constraint ( 14) is removed, model ( 23) can be converted into a multiple traveling salesman problem with time window (m-TSPTW).

Example Setting
In this example, the logistics network consisted of a harvest location (0) and 15 consumers (1-15).The harvest location had six homogeneous vehicles with the average speed of 60 km/h.The harvest location provided one kind of perishable agri-product.The unit harvesting time was 0.007 h for a unit perishable agri-product.The quality decay rate was 0.5 unit/h.The time window of harvest location was [8:00, 17:00] (For simplicity without generality, 8 : 00 was assumed as 0, and 17:00 was assumed as 9 during the calculation in the example.).Table 1 gives the delivery time between each node located in different geographical locations, in which M denotes an infinite number, which refers to the node pair not being linked geographically.Figure 3 demonstrates the spatial location of each node in the coordinate according to the data given in Solomon's problem [40].The coordinate of the harvest location was (17.5, 17.5).Table 2 gives the spatial location, demand, and time windows of each consumer.

Numerical Results
The numerical results are reported in Table 3.It can be seen that the optimal route of vehicle 1 starts from harvest location 0, and serves the consumers with the sequence 10 → 5 → 9 → 2 → 6, and then returns to harvest location 0. The optimal starting harvest time for the distribution assignments by vehicle 1 is 8: 50 a.m.Similarly, the optimal route of vehicle 2 is 0 → 8 → 14 → 12 → 11 → 4 → 0 , with the optimal starting harvest time for the distribution assignments by vehicle 2 being 9: 44 a.m.For vehicle 3, the optimal route is 0 → 7 → 15 → 1 → 3 → 13 → 0, with the optimal starting harvest time 8: 14 a.m.According to the results given in Table 3, Figure 4 depicts the optimal vehicle route of vehicle 1, vehicle 2, and vehicle 3. Note that, the optimal vehicle route of vehicle 1 is 0 → 10 → 5 → 9 → 2 → 6 → 0. It can obviously be seen that the total route length of the optimal route is much longer than the vehicle routing sequence 0 → 10 → 5 → 2 → 9 → 6 → 0. That is because the time window of consumer 9 is earlier than consumer 2, which constrains vehicle 1 to serve consumer 9 first rather than according to the shortest path.Similarly, the optimal vehicle route of vehicle 2 presents the same phenomena.Therefore, the results given in Table 3 and Figure 4 show that the optimal routes are not the shortest paths due to the impact of time window constraint in the proposed model.

Numerical Results
The numerical results are reported in Table 3.It can be seen that the optimal route of vehicle 1 starts from harvest location 0, and serves the consumers with the sequence 10 → 5 → 9 → 2 → 6 , and then returns to harvest location 0. The optimal starting harvest time for the distribution assignments by vehicle 1 is 8:50 a.m.Similarly, the optimal route of vehicle 2 is 0 → 8 → 14 → 12 → 11 → 4 → 0 , with the optimal starting harvest time for the distribution assignments by vehicle 2 being 9:44 a.m.For vehicle 3, the optimal route is 0 → 7 → 15 → 1 → 3 → 13 → 0 , with the optimal starting harvest time 8:14 a.m.According to the results given in Table 3, Figure 4 depicts the optimal vehicle route of vehicle 1, vehicle 2, and vehicle 3. Note that, the optimal vehicle route of vehicle 1 is 0 → 10 → 5 → 9 → 2 → 6 → 0 .It can obviously be seen that the total route length of the optimal route is much longer than the vehicle routing sequence 0 → 10 → 5 → 2 → 9 → 6 → 0. That is because the time window of consumer 9 is earlier than consumer 2, which constrains vehicle 1 to serve consumer 9 first rather than according to the shortest path.Similarly, the optimal vehicle route of vehicle 2 presents the same phenomena.Therefore, the results given in Table 3 and Figure 4 show that the optimal routes are not the shortest paths due to the impact of time window constraint in the proposed model.
6 → 0. It can obviously be seen that the total route length of the optimal route is much longer than the vehicle routing sequence 0 → 10 → 5 → 2 → 9 → 6 → 0. That is because the time window of consumer 9 is earlier than consumer 2, which constrains vehicle 1 to serve consumer 9 first rather than according to the shortest path.Similarly, the optimal vehicle route of vehicle 2 presents the same phenomena.Therefore, the results given in Table 3 and Figure 4 show that the optimal routes are not the shortest paths due to the impact of time window constraint in the proposed model.

Sensitivity Analysis
We further carried out sensitivity analysis in terms of different quality decay rate, different time windows, and different number of vehicles.

Sensitivity Analysis on Quality Decay Rate
We varied the quality decay rate from 0.1 to 0.9.The numerical results are reported in Table 4 and the objective values are given after normalization. Figure 5 depicts the relationship between the quality decay rate and the model objective value.It can be seen that the objective value is increasing with the growing of the quality decay rate.In other words, if the perishable agri-products has a higher quality decay rate, the harvest location will adopt more vehicles to distribute the perishable agri-products to consumers, so as to maintain the freshness and reduce the deterioration loss.

Sensitivity Analysis
We further carried out sensitivity analysis in terms of different quality decay rate, different time windows, and different number of vehicles.

Sensitivity Analysis on Quality Decay Rate
We varied the quality decay rate from 0.1 to 0.9.The numerical results are reported in Table 4 and the objective values are given after normalization. Figure 5 depicts the relationship between the quality decay rate and the model objective value.It can be seen that the objective value is increasing with the growing of the quality decay rate.In other words, if the perishable agri-products has a higher quality decay rate, the harvest location will adopt more vehicles to distribute the perishable agriproducts to consumers, so as to maintain the freshness and reduce the deterioration loss.

Sensitivity Analysis on Time Windows
We changed the number of vehicles from three to nine and relaxed the time window constraint under four cases, denoted as TW 100% , TW 80% , TW 60% , and TW 40% .TW 100% means that all consumers have time window requirements; TW 80% means that 80% of consumers have time window requirements; TW 60% means that 60% of consumers have time window requirements; and TW 40% means that 40% of consumers have time window requirements.Figure 6 shows that with the increasing number of vehicles, the objective value improves under four cases of TW requirements and is almost approaching the same value if the vehicle number equals nine.We conclude that the objective value will improve rapidly when the vehicle number increases more because of the large proportion of fixed costs in objective value.

Sensitivity Analysis on Time Windows
We changed the number of vehicles from three to nine and relaxed the time window constraint under four cases, denoted as TW (100%), TW (80%), TW (60%), and TW (40%).TW (100%) means that all consumers have time window requirements; TW (80%) means that 80% of consumers have time window requirements; TW (60%) means that 60% of consumers have time window requirements; and TW (40%) means that 40% of consumers have time window requirements.Figure 6 shows that with the increasing number of vehicles, the objective value improves under four cases of TW requirements and is almost approaching the same value if the vehicle number equals nine.We conclude that the objective value will improve rapidly when the vehicle number increases more because of the large proportion of fixed costs in objective value.

Sensitivity Analysis on Vehicles
We changed the number of the vehicles from three to nine and calculated the average loading ratio under the time window constraint of each consumer.The average loading ratio refers to the average value of all vehicles' loading ratio, while the loading ratio equals the total vehicle distributed amount divided by its capacity. Figure 7 depicts the variation trend between the average loading ratio and the number of vehicles.Obviously, with an increasing number of vehicles, the average loading ratio goes down.According to the above discussions, we can summarize two management insights from this study.On one hand, increasing the number of vehicles could reduce the quality decay of perishable agri-products, and at the same time, it could also lower the average loading ratio and increase the total cost.On the other hand, in order to meet the fresh quality requirement of consumers, the harvest location could increase the number of vehicles to prevent quality deterioration of perishable agriproducts.

Conclusions
In this study, we have proposed mixed integer nonlinear programming to model the integrated problem of perishable agri-products harvest and distribution scheduling with time windows, with the objectives of minimizing the average quality decay cost and the distribution cost.In the proposed model, we initially proposed an average quality decay function to capture the quality decay of perishable agri-products during the harvest and distribution phases.Therefore, a big M method was

Sensitivity Analysis on Vehicles
We changed the number of the vehicles from three to nine and calculated the average loading ratio under the time window constraint of each consumer.The average loading ratio refers to the average value of all vehicles' loading ratio, while the loading ratio equals the total vehicle distributed amount divided by its capacity. Figure 7 depicts the variation trend between the average loading ratio and the number of vehicles.Obviously, with an increasing number of vehicles, the average loading ratio goes down.

Sensitivity Analysis on Vehicles
We changed the number of the vehicles from three to nine and calculated the average loading ratio under the time window constraint of each consumer.The average loading ratio refers to the average value of all vehicles' loading ratio, while the loading ratio equals the total vehicle distributed amount divided by its capacity. Figure 7 depicts the variation trend between the average loading ratio and the number of vehicles.Obviously, with an increasing number of vehicles, the average loading ratio goes down.According to the above discussions, we can summarize two management insights from this study.On one hand, increasing the number of vehicles could reduce the quality decay of perishable agri-products, and at the same time, it could also lower the average loading ratio and increase the total cost.On the other hand, in order to meet the fresh quality requirement of consumers, the harvest location could increase the number of vehicles to prevent quality deterioration of perishable agriproducts.

Conclusions
In this study, we have proposed mixed integer nonlinear programming to model the integrated problem of perishable agri-products harvest and distribution scheduling with time windows, with the objectives of minimizing the average quality decay cost and the distribution cost.In the proposed model, we initially proposed an average quality decay function to capture the quality decay of perishable agri-products during the harvest and distribution phases.Therefore, a big M method was adopted to convert the nonlinear objective function into a linear form.In addition, numerical According to the above discussions, we can summarize two management insights from this study.On one hand, increasing the number of vehicles could reduce the quality decay of perishable agri-products, and at the same time, it could also lower the average loading ratio and increase the total cost.On the other hand, in order to meet the fresh quality requirement of consumers, the harvest location could increase the number of vehicles to prevent quality deterioration of perishable agri-products.

Conclusions
In this study, we have proposed mixed integer nonlinear programming to model the integrated problem of perishable agri-products harvest and distribution scheduling with time windows, with the objectives of minimizing the average quality decay cost and the distribution cost.In the proposed model, we initially proposed an average quality decay function to capture the quality decay of perishable agri-products during the harvest and distribution phases.Therefore, a big M method was adopted to convert the nonlinear objective function into a linear form.In addition, numerical experiments have been conducted to verify the efficiency and effectiveness of the proposed model, and the numerical results and sensitivity analysis have also been reported.Finally, we discussed how the results could serve as a guide for operations management of agri-products in practice.
There are specific areas in this study that can be researched further.This study only considered a single harvest location with a single type of perishable agri-product, hence multiple harvest locations and multiple types of perishable agri-products could be considered as a prospect.Another prospect would consider the penalty cost for time windows in the objective function.Finally, future studies could incorporate more uncertain factors into the integrated harvest and distribution scheduling, such as the uncertainties of consumer demand and time windows, and the uncertainties of temperature and weather condition.We leave these problems for the future work.

Figure 1 .
Figure 1.Schematic graph of harvest and distribution scheduling.

Figure 2 .
Figure 2. Quantity variation of consumer i during the harvest time period.

Figure 1 .
Figure 1.Schematic graph of harvest and distribution scheduling.

Figure 1 .
Figure 1.Schematic graph of harvest and distribution scheduling.

Figure 2 .
Figure 2. Quantity variation of consumer i during the harvest time period.

Figure 2 .
Figure 2. Quantity variation of consumer i during the harvest time period.

Figure 4 .
Figure 4. Optimal route of three vehicles.Figure 4. Optimal route of three vehicles.

Figure 4 .
Figure 4. Optimal route of three vehicles.Figure 4. Optimal route of three vehicles.

Figure 5 .
Figure 5. Relationship between the quality decay rate and the objective value.

Figure 5 .
Figure 5. Relationship between the quality decay rate and the objective value.

Sustainability 2018 , 13 Figure 6 .
Figure 6.Objective value under different number of vehicles and time windows.

Figure 7 .
Figure 7.The variation trend between average loading ratio and number of vehicles.

Figure 6 .
Figure 6.Objective value under different number of vehicles and time windows.

Sustainability 2018 , 13 Figure 6 .
Figure 6.Objective value under different number of vehicles and time windows.

Figure 7 .
Figure 7.The variation trend between average loading ratio and number of vehicles.

Figure 7 .
Figure 7.The variation trend between average loading ratio and number of vehicles.

Table 1 .
Delivery time for the node pairs (h).

Table 2 .
Demand requirement and time windows of each node.
Note: Consumer number is abbreviated as Con.No.

Table 3 .
Distribution routings of each vehicle.

Table 3 .
Distribution routings of each vehicle.

Table 4 .
Objective value for different quality decay rate.

Table 4 .
Objective value for different quality decay rate.