TY - JOUR AU1 - Guitián, de Frutos, Roque M AU2 - Casas-Méndez,, Balbina AB - Abstract This study considers the problem of optimizing the routes of vehicles used in an agricultural cooperative that distributes animal feed to its customers. Different peculiarities distinguish our problem from the general Vehicle Routing Problem because there exists a type of time window constraint, truck capacity is limited, trucks are compartmentalized, there are access restrictions for trucks to some farms and our problem has the main objective of maximizing the amount of feed distributed every day and an auxiliary criterion of minimizing the transport costs. We present a mathematical programming formulation of the model that represents the logistics management problem. As solving the exact model is computationally expensive, we opt for a hybrid heuristic approach that first uses an insertion heuristic algorithm to obtain an initial solution and then improvement processes are performed using the so-called simulated annealing metaheuristic. We have built a graphical interface that allows easy use of the system and, in particular, efficient post-optimality analysis. Moreover, this interface can easily interface with other management tools. The utility of our model is shown with real and simulated data sets. 1. Introduction The logistics of transport is a subject of increasing importance as societies advance their industrial development. The use of computational methods to calculate optimal or near-optimal routes can result in remarkable savings. This is why public and private companies are interested in developing flexible tools that automate the design of routes with additional objectives such as minimizing costs. Transportation management often follows a pattern in which a centralized fleet distributes goods from a central warehouse to a set of receivers that are positioned in different locations. As the number of vehicles, deliveries, goods or receivers grows, so do the difficulties in arranging a plan that minimizes costs. This work is motivated by a real problem from a company focused on the manufacture and distribution of feed. AIRA is the name of an agricultural cooperative (formed 12 years ago) in Taboada, a Galician village located halfway between the cities of Lugo and Orense, in northwest Spain. It produces four different types of feed that are distributed to its customers. The company has approximately 1500 farmers (clients) scattered over 60 adjacent municipalities, covering a wide area, as shown in Fig. 1. Fig. 1. Open in new tabDownload slide Geographic area of influence of the Spanish agricultural cooperative AIRA. Fig. 1. Open in new tabDownload slide Geographic area of influence of the Spanish agricultural cooperative AIRA. Farmers place orders for different types of feed, which may be urgent or have a specified delivery time. For urgent orders, the period is reduced to only 1 day and such orders cost a higher price. Farmers usually issue orders once or twice a month and, under normal conditions, only require one type of feed. The distribution of feed is performed by a fleet of trucks and each driver is paid according to the distance travelled and the load. The trucks are different sizes and legal requirements limit the transportable mass and distance travelled each day. The vehicles are compartmentalized into several hoppers (three to five, generally), with each having a different capacity between 1700 and 4500 kg. Moreover, for technical reasons, each bin can contain only one type of food and cannot be shared between farmers. Moreover, given the size and characteristics of each truck, some trucks cannot travel on certain pathways or access certain farms, so there is a list of customers who cannot be served by each vehicle. The amount of feed produced in the past year was 146000. As a sample of the type of orders usually received, in Fig. 2 one can see the volume of 36 orders (urgent or not) planned on a certain day. As can be seen, there are differences among the orders from that day and among the orders on different days. Briefly, the average order of a partner in a workday is between 2700 and 8700 kg, the average annual order per member is between 46000 and 148000 kg and the average number of annual orders of a partner is around 17. However, all of these parameters vary each year, with the disappearance of some livestock farms, mergers among cooperatives, the arrival of new partners or other eventualities. Fig. 2. Open in new tabDownload slide Volume (kg) of the 36 orders placed by the members of the cooperative on a certain day. Fig. 2. Open in new tabDownload slide Volume (kg) of the 36 orders placed by the members of the cooperative on a certain day. The purpose of this study is to provide a tool for the cooperative that automatically proposes routes for each truck while satisfying the constraints and that maximizes the amount of feed distributed daily while minimizing transport costs. These transport costs are broken down into a fixed component for each truck, plus the sum of an amount proportional to the number of deliveries on each route and an amount that is proportional to the distance travelled and amount of feed carried on each route. This means that one of the objective functions will be non-linear in our model. It is convenient to indicate that on any day there may be urgent orders, the trucks may experience mechanical failure and the forecasts about the level of feed on farms may not be reliable or the partners may communicate their needs at short notice. All of these factors suggests that, at the moment, the designed tool is especially useful for planning the routes for 2 or 3 days ahead at most. A research group of Agricultural Engineering at the University of Santiago de Compostela (USC) has designed a comprehensive global positioning system (GPS) that can monitor various truck routes (cf. Amiama-Ares et al., 2014). However, cooperative managers hope to improve the system in the coming years. Hence, this routing work is a contribution to the automatization of the distribution process of the company. The cooperative also plans to equip the silos in which feed is stored with sensors, making it possible to remotely collect data on available volume and daily consumption. Thus, the system may be able to indicate the silos that require restocking, the quantity of feed required and the deadline. Finally, automation of truck hopper opening in the silo charging process is also planned. The currently used system (GPS) provides all the geographical information needed to solve the problem. As for the data used, we know the capacity of each bin, the demands of different customers and the urgencies of the orders. We also know whether a truck can access each farm and its maximum load. Moreover, we use the information regarding the price of each download and the variable cost of transporting a certain amount of food over a given distance. 2. Theoretical framework and literature review The problem in question is a modification of the well-known Vehicle Routing Problem (VRP), which is a combinatorial optimization problem. This kind of problem involves trying to minimize or maximize an objective function of several variables defined on a discrete set. The interest in this area of research is both theoretical and applied. Most VRPs belong to the complexity class NP (nondeterministic polynomial time)-complete. The academic motivation to solve these problems is that it is not possible to construct a polynomial time algorithm that solves any instance of the problem. The VRPs were introduced by Dantzig & Ramser (1959). In our problem, we have an heterogeneous fleet of trucks that can travel several routes, in each working day, to meet the demands of the cooperative members. From a theoretical point of view, different peculiarities distinguish our problem from the general VRP. First, there is a kind of time window constraint, since some partners need to receive the food order urgently. In this sense, our problem is a generalization of the model of Golden & Assad (1986). Second, it is noteworthy that the truck capacity is limited as described in the line of the work of Labbé et al. (1991). Third, the trucks are compartmentalized. This allows partner demand to be met by dividing the order among several compartments. For this reason, as in the model of a VRP with split delivery (Dror & Trudeau, 1989) we included continuous variables in our model indicating the proportion of the request of a partner served in a compartment (of a truck, in one of its routes). The work of Derigs et al. (2011) is another example where routing problems are considered in which vehicles are divided into compartments. It should be mentioned that in our case, for technical reasons, a compartment cannot be used to meet the demand of two different partners and can only be occupied by one type of feed. Fourth, there are restrictions of accessibility to trucks on certain farms. This restriction is considered in Chao (2002), Hoff (2006) and Derigs et al. (2013), among others. Indeed, these authors considered the so-called truck and trailer routing problems, where a fleet of vehicles consisting of trucks and trailers attends a set of customers. When a trailer cannot access a certain customer, this is attended only by the truck. In our case, we do not work with trailers, so if a truck cannot access a client, it must necessarily be served by a different truck. Fifth, and finally, our problem has a main objective of maximizing the number of members served each day, but also an auxiliary criterion that is given by the desire to minimize costs. Tavakkoli-Moghaddam et al. (2007) also considered an additional objective to minimize costs, which in that case was to maximize the utilization capacity of the trucks. The simultaneous combination of these five ingredients makes our problem singular. Heuristic algorithms have been effective tools for solving NP-hard problems (cf. Osman & Laporte, 1996). Dror & Trudeau (1989) proposed a local search in two stages to solve the VRP with split delivery. Mullaseril et al. (1997) adapted this algorithm to solve a problem of feed distribution on a cattle ranch in Arizona, where there are also time window restrictions. Bouzaiene-Ayari et al. (1993) proposed a modification of the algorithm of Clarke & Wright (1964) to solve the VRP with stochastic demands and split delivery. Perboli et al. (2017) used a heuristic method to solve the so-called multipath travelling salesman problem with stochastic travel costs, specifically designed for smart cities and city logistics applications. Ho & Haugland (2004) designed a tabu search algorithm for the VRP with split delivery and time window restrictions. Carpente et al. (2010) used a tabu search algorithm to solve the problem of the harvesters of an agricultural cooperative that must cross the different smallholdings of the partners according to temporary requests and so that the smallholdings of the same partner must be treated in block for technical reasons. Ruiz et al. (2004) were also inspired by the real problem of another Spanish cooperative that manufactures animal feed. Their approach generates all the feasible routes by means of an implicit enumeration algorithm as a first step and then used an integer programming model that selects the optimal routes. They reduced the initial complexity of the problem by clustering the customers that they have to serve. Ambrosino & Sciomachen (2007) considered a variant of the VRP with capacities and split delivery using heuristic methods that solve the problem in two stages. They started by making clusters with all customers and they continued with an improvement algorithm through local search. It was applied to an Italian company that is facing a problem of food distribution. A case of food distribution in Ukraine was studied by Kuznietsov et al. (2016). The problem was heuristically broken into subproblems, based on a partition of the geographical area where there are customers, to optimize routes and reduce the number of vehicles used. Chen et al. (2007) combined a mixed integer program and an algorithm that starts with an initial solution for the VRP with split delivery, also based on the algorithm of Clarke & Wright (1964). Brito et al. (2016) analysed a variant of the VRP with time windows where the company has its own fleet of vehicles and also hires the service of other companies that provide complementary resources for specific service needs. They used a variable neighbourhood search algorithm to solve this problem. Other related problems are dynamic combinatorial optimization problems that are those in which decision variables, objectives and/or constraints can change over time. Yang et al. (2013) considered bio-inspired metaheuristics in the context of such problems. Simulated annealing (cf. Kirkpatrick et al., 1983; Nunes De Castro, 2006) is a heuristic based on local search that deviates from local optima accepting in its iterations worse solutions with a small probability. It has a beautiful physical interpretation and its simplicity makes it competitive computationally with other types of metaheuristic algorithms. In addition, it has been successfully used in many applications. Bent & Van Hentenryck (2004) used a local search for the VRP with time windows that combined simulated annealing and large neighbourhood search to first minimize the number of vehicles and then minimize travel costs. They tested the algorithm using benchmarks from the literature. Tavakkoli-Moghaddam et al. (2007) also used a simulated annealing algorithm to solve a problem with heterogeneous vehicles and they tested the algorithm on randomly generated instances. Lin et al. (2011) used a simulated annealing algorithm for solving a routing problem of trucks and trailers with time windows and they generated tests for a computational study and tested the algorithm. 3. Methodological approach As a consequence of the above considerations, we propose a formulation of this problem using a nonlinear model, which we write in the language A Mathematical Programming Language (AMPL) (cf. http://ampl.com/). Models prepared in AMPL can be run on the NEOS server (cf. https://neos-server.org/neos/) or using a free or commercial license. For small instances, we contrast and solve the model using a lexicographic procedure (cf. Hwang & Masud, 1979, for more details about lexicographic optimization methods in programming problems with several objectives). In this way, we confront the task of adapting more classical vehicle routing models to the specific case of AIRA. Then, we use an algorithm that provides an ‘exact solution’. The exact resolution of small examples with two trucks and six clients, by means of the NEOS optimization solver, requires more than 3 h of computing time. Therefore, we looked for an alternative procedure for practical use. In this way, we opted for a hybrid heuristic approach that first uses an insertion heuristic algorithm to obtain an initial solution (cf. Clarke & Wright, 1964; Mole & Jameson, 1976). Then, improvement processes are performed on the initial solution using a simulated annealing metaheuristic. With this procedure, if we consider medium or large instances, involving more than 60 orders and 6 trucks, we obtain an acceptable solution in less than 10 min. As a complement, we have created a graphical user interface in JAVA language that allows data to be introduced, the results to be viewed in detail and the routes obtained to be drawn on a map, and it provides a general user-friendly interface for the technicians of the company. Finally, by means of real and simulated data, different types of tests have been performed that guarantee the adequacy of the tool created, both from the point of view of its efficiency and the quality of the solutions it allows to be obtained. 4. Formulation of the problem In this section, we present a mathematical programming formulation of the model that represents the logistic management problem of AIRA. 4.1 Notation and decision variables The main elements of our formulation are as follows. – |$\mathscr{N}$| represents the set of stockbreeders or members of the cooperative. – |$\overline{\mathscr{N}}=\mathscr{N}\cup \{0\}$| is the set of stockbreeders and the cooperative itself. – |$\mathscr{T}$| is the set of available trucks. – For each |$t\in \mathscr{T}$|⁠, |$\mathscr{R}_t$| is the set of possible routes of truck |$t$| that are considered in the planning. – For each |$t\in \mathscr{T}$|⁠, |$\mathscr{H}_t$| is the set of hoppers of truck |$t$|⁠. – For each |$t\in \mathscr{T}$|⁠, let |$C^t=(c_h^t)_{h\in \mathscr{H}_t}$| be the vector of capacities of truck |$t$|⁠, where |$c_h^t$| is the capacity of hopper |$h\in \mathscr{H}_t$|⁠. – |$\mathscr{F}$| is the set of livestock feeds available. – |$O=(o_{n,f})_{n\in \mathscr{N},f\in \mathscr{F}}$| is the matrix of orders, where |$o_{n,f}$| represents the order placed by |$n\in \mathscr{N}$| for feed type |$f\in \mathscr{F}$|⁠. – |$D=(d_{n_1,n_2})_{(n_1,n_2)\in \overline{\mathscr{N}}\times \overline{\mathscr{N}}}$| is the matrix of distances in minutes, i.e. |$d_{n_1,n_2}$| is the time a truck needs to go from |$n_1$| to |$n_2$|⁠. – |$E=(e_{n_1,n_2})_{(n_1,n_2)\in \overline{\mathscr{N}}\times \overline{\mathscr{N}}}$| is the matrix of distances, i.e. |$e_{n_1,n_2}$| is the distance in kilometres from |$n_1$| to |$n_2$|⁠. – |$I=(i_{t,n})_{t\in \mathscr{T},n\in \mathscr{N}}$| is the matrix that describes which stockbreeders each truck can visit, where |$i_{t,n}$| equals |$1$| if truck |$t$| can visit client |$n$| and is |$0$| otherwise. – |$L=(l_t)_{t\in \mathscr{T}}$| is the maximum allowable load that truck |$t$| is legally authorized to carry. – |$U=(u_n)_{n\in \mathscr{N}}$| is the urgency vector where for each |$n\in \mathscr{N}$|⁠, |$u_n$| denotes the maximum number of days that stockbreeder |$n$| can wait to be served. – |$M$| is the estimated maximum distance that a truck can cover in a single day. – |$\alpha $| is the fixed cost for each gramme carried. We include this parameter in the notation for completeness. However, we will not find it in the formulation of restrictions or objectives because it does not affect them. As we want, at one time or another, to attend to all orders, this generates a fixed cost of |$\alpha $| multiplied by the total grammes demanded, regardless of the method of designing the routes. – |$\beta $| is the fixed cost of each unloading of a truck. – |$\gamma (a)$| is the variable cost per gramme carried and kilometre travelled, where |$a$| is the distance covered. Two sets of decision variables are considered. The first is composed of the variables, |$x_{n_1,n_2}^{t,r}$|⁠, given by $$x_{n_1,n_2}^{t,r}= \begin{cases} 1 \quad & \text{if truck} \ \textit{t}, \ \text{on its route} \textit{r}, \text{travels from} \ n_1 \mbox{ to } n_2,\\ 0 \quad & \mbox{otherwise,}\end{cases}$$ for each |$t\in \mathscr{T}$|⁠, |$r\in \mathscr{R}_t$| and |$n_1,n_2\in \overline{\mathscr{N}}$|⁠. The second set of variables, |$y_{n,f,h}^{t,r}\in [0,1]$|⁠, denotes for each truck |$t\in \mathscr{T}$| the percentage of the hopper |$h\in \mathscr{H}_t$| that is being used to carry feed |$f\in \mathscr{F}$| on the route |$r\in \mathscr{R}_t$| to serve the stockbreeder |$n\in \mathscr{N}$|⁠. The first set of variables determines the routes that each truck must follow in its working day and the second set of variables represents the way in which the trucks must be loaded. The need for a continuous set of variables is motivated by the form of the objective functions that the mathematical programming problem must optimize, which are introduced in Section 4.3. 4.2 Constraints The relationships that describe the real-world model are translated into our formulation via mathematical constraints. As is usual in models of a certain complexity, we include restrictions comprehensively. Even though some may seem trivial or redundant, we add them to the model to reinforce the formulation with the aim of achieving better computational results. The constraints are grouped into five sets. Constraints describing the routes of trucks, Constraints that link the routes of trucks and their loading processes, Constraints modelling the requests of the stockbreeders on the urgency of the orders (a type of time-window constraint), Constraints modelling the capacities of trucks and limits on the hours and kilometres in the working day (capacity constraints) and Constraints modelling technical restrictions on the procedure for loading the trucks. We now describe each of the constraints in detail. Constraints describing the routes of trucks. No truck is allowed to go from a stockbreeder and back to the same stockbreeder: $$\begin{equation*}x_{n,n}^{t,r}=0 \quad \forall\ t\in\mathscr{T},\quad \forall\ r\in \mathscr{R}_t,\ \quad \forall\ n\in\overline{\mathscr{N}}.\end{equation*}$$ Each route of each truck departs from the cooperative headquarters at most once: $$\begin{equation*}\sum_{n\in\mathscr{N}}x_{0,n}^{t,r}\leqslant 1 \quad \forall\ t\in\mathscr{T},\ \quad\forall\ r\in \mathscr{R}_t.\end{equation*}$$ Each planned route of each truck has to start from the cooperative headquarters: $$\begin{equation*}x_{n_1,n_2}^{t,r}\! \leqslant \!\sum_{n\in\mathscr{N}}x_{0,n}^{t,r} \quad \forall\ t\in\mathscr{T},\quad \forall\ r\in \mathscr{R}_t,\quad \forall\ n_1,n_2\in\mathscr{N}.\end{equation*}$$ Flow conservation constraint. On each route, if a truck arrives at a stockbreeder, it also has to leave from the stockbreeder: $$\begin{equation*}\sum_{n\in\overline{\mathscr{N}}}x_{n,n_1}^{t,r}=\sum_{n\in\overline{\mathscr{N}}}x_{n_1,n}^{t,r} \quad \forall\ t\in\mathscr{T}, \quad \forall\ r\in \mathscr{R}_t,\ \quad \forall\ n_1\in\overline{\mathscr{N}}.\end{equation*}$$ Subtour elimination on each route of each truck: $$\begin{equation*}\sum_{n_1,n_2\in S}x_{n_1,n_2}^{t,r}\leqslant\ |S|-1 \quad \forall\ t\in \mathscr{T},\quad \forall\ r\in\mathscr{R}_t,\quad \forall\ S\subseteq\mathscr{N} : |S|>1.\end{equation*}$$ Constraints that link the routes of trucks and their loading. If a truck is carrying feed for a stockbreeder, we ensure that the truck visits the stockbreeder: $$\begin{equation*}\frac{1}{\mid \mathscr{H}_t \mid}\sum_{f\in \mathscr{F}}\sum_{h\in \mathscr{H}_t} y_{n_1,f,h}^{t,r} \leqslant\ \sum_{n\in\overline{\mathscr{N}}}x_{n,n_1}^{t,r}\quad \forall\ t\in\mathscr{T},\quad \forall\ r\in \mathscr{R}_t,\quad \forall\ n_1\in \mathscr{N}.\end{equation*}$$ If a truck visits a stockbreeder on a route, it is indeed carrying feed for him: $$\begin{equation*}\sum_{n\in\overline{\mathscr{N}}}x_{n,n_1}^{t,r}\leqslant\ \sum_{f\in \mathscr{F}}\sum_{h\in \mathscr{H}_t} c_h^t \cdot y_{n_1,f,h}^{t,r}\quad \forall\ t\in\mathscr{T},\quad \forall\ r\in \mathscr{R}_t,\quad \forall\ n_1\in \mathscr{N}.\end{equation*}$$ Constraints modelling urgent stockbreeder orders (a type of time-window constraint). Stockbreeders with urgent orders are served their entire demanded orders: $$\begin{equation*}\sum_{t\in\mathscr{T}}\sum_{r\in\mathscr{R}_t}\sum_{h\in \mathscr{H}_t}c_h^t \cdot y_{n,f,h}^{t,r}=o_{n,f}\quad \forall\ f\in \mathscr{F},\quad \forall\ n\in\mathscr{N} : u_n=0.\end{equation*}$$ No stockbreeder is served more than what they actually ordered: $$\begin{equation*}\sum_{t\in\mathscr{T}}\sum_{r\in\mathscr{R}_t}\sum_{h\in \mathscr{H}_t}c_h^t \cdot y_{n,f,h}^{t,r}\leqslant\ o_{n,f}\quad \forall\ f\in \mathscr{F},\quad \forall\ n\in\mathscr{N} : u_n>0.\end{equation*}$$ Constraints modelling the capacities of trucks and limits on hours and kilometres in the working day (capacity constraints). Impossible visits are not allowed: $$\begin{equation*}{\textrm{In maximization problem,}}\quad \sum_{n\in\overline{\mathscr{N}}}x_{n,n_1}^{t,r}\leqslant i_{t,n_1}\quad \forall\ t\in\mathscr{T},\quad \forall\ r\in\mathscr{R}_t,\quad \forall\ n_1\in\mathscr{N}.\end{equation*}$$ $$\begin{equation*}{\textrm{In minimization problem,}}\quad \sum_{f\in \mathscr{F}}y_{n_1,f,h}^{t,r}\leqslant i_{t,n_1}\quad \forall\ t\in\mathscr{T},\quad \forall\ r\in\mathscr{R}_t,\quad \forall\ h\in \mathscr{H}_t,\quad\!\! \forall\ n_1\in\mathscr{N}.\end{equation*}$$ The working day of each truck lasts a maximum of 9 h (540 min). We do not consider unloading time: $$\begin{equation*}\sum_{r\in\mathscr{R}_t}\sum_{n_1\in\overline{\mathscr{N}}}\sum_{n_2\in\overline{\mathscr{N}}} d_{n_1,n_2} \cdot x_{n_1,n_2}^{t,r}\leqslant 540\quad \forall\ t\in\mathscr{T}.\end{equation*}$$ A truck cannot travel more than a certain number of kilometres in a day: $$\begin{equation*}\sum_{r\in\mathscr{R}_t}\sum_{n_1\in\overline{\mathscr{N}}}\sum_{n_2\in\overline{\mathscr{N}}} e_{n_1,n_2} \cdot x_{n_1,n_2}^{t,r}\leqslant M\quad \forall\ t\in\mathscr{T}.\end{equation*}$$ No truck carries more than the maximum load that it is legally authorized to carry: $$\begin{equation*}\sum_{n\in\mathscr{N}}\sum_{f\in\mathscr{F}}\sum_{h\in \mathscr{H}_t}c_h^t \cdot y_{n,f,h}^{t,r}\leqslant l_t\quad \forall\ t\in \mathscr{T},\quad \forall\ r\in\mathscr{R}_t.\end{equation*}$$ Constraints modelling technical restrictions on the loading procedure of the trucks. The hopper of each truck cannot be used to serve more than one stockbreeder: $$ y_{n_1,f_1,h}^{t,r}+y_{n_2,f_2,h}^{t,r}\leqslant\ \max\left\{y_{n_1,f_1,h}^{t,r},y_{n_2,f_2,h}^{t,r}\right\}\quad \forall\ t\in \mathscr{T},\quad \forall\ r\in\mathscr{R}_t,\quad \forall\ h\in \mathscr{H}_t,\ $$ $$ \forall\ n_1,n_2\in \mathscr{N},\quad \forall\ f_1,f_2\in\mathscr{F} : n_1\neq n_2\ \textrm{and} f_1\neq f_2; $$ $$\begin{equation*}y_{n_1,f,h}^{t,r}+y_{n_2,f,h}^{t,r}\leqslant\ \max\left\{y_{n_1,f,h}^{t,r},y_{n_2,f,h}^{t,r}\right\}\quad \begin{array}{c} \forall\ t\in \mathscr{T},\ \forall\ r\in\mathscr{R}_t,\ \forall\ h\in \mathscr{H}_t, \\ \forall\ f\in\mathscr{F},\ \forall\ n_1,n_2\in \mathscr{N} : n_1\neq n_2\ .\end{array}\end{equation*}$$ The hopper of each truck cannot be loaded with more than one type of feed: $$\begin{equation*}y_{n,f_1,h}^{t,r}+y_{n,f_2,h}^{t,r}\leqslant\ \max\left\{y_{n,f_1,h}^{t,r},y_{n,f_2,h}^{t,r}\right\}\quad \begin{array}{c} \forall\ t\in \mathscr{T},\ \forall\ r\in\mathscr{R}_t,\ \forall\ h\in \mathscr{H}_t, \\ \forall\ f_1,f_2\in\mathscr{F},\ \forall\ n\in \mathscr{N} : f_1\neq f_2 .\end{array}\end{equation*}$$ 4.3 Objective function As mentioned above, we face a programming model with a main objective and a second auxiliary objective. First, we maximize the amount of feed delivered in the current working day. In other words, we first try to serve as many stockbreeders with non-urgent orders as possible. Recall that the first constraint in group 3 implies that stockbreeders with urgent orders are served within the requested times. As a second objective, the total transportation cost assumed by the cooperative is minimized. Thus, $$\begin{align*}&\mbox{first,} \qquad \max\left\{\sum_{\substack{t\in\mathscr{T}\\r\in \mathscr{R}_t\\h\in\mathscr{H}_t}}\sum_{\substack{f\in\mathscr{F}\\n\in\mathscr{N}}}c_h^t\cdot y_{n,f,h}^{t,r}\right\}, \qquad \mbox{and then}\\ &\min \left\{\sum_{\substack{t\in\mathscr{T}\\r\in \mathscr{R}_t}}\left[\vphantom{\sum_{\substack{h\in\mathscr{H}_t\\f\in\mathscr{F},n\in\mathscr{N}}}}\beta\sum_{n_1,n_2\in\overline{\mathscr{N}}}x_{n_1,n_2}^{t,r} \right.\right.\left.\left. +\ \gamma \left(\sum_{n_1,n_2\in\overline{\mathscr{N}}}e_{n_1,n_2}\cdot x_{n_1,n_2}^{t,r}\right)\!\sum_{n_1,n_2\in\overline{\mathscr{N}}}\!e_{n_1,n_2}\cdot x_{n_1,n_2}^{t,r}\!\sum_{\substack{h\in\mathscr{H}_t\\f\in\mathscr{F},n\in\mathscr{N}}}\!c_h^t\cdot y_{n,f,h}^{t,r}\right]\right\}.\end{align*}$$ It should be noted that the present model is valid for designing the 1-day route. Consequently, to plan several days, it is necessary to execute it several times, in order, so that for each day, we will consider the corresponding urgency of the orders and the previous results. There could be an alternative model that takes into account, when maximizing the load, the urgency of all orders. However, according to the experience of the company, this is not necessary, and therefore, we discard this idea, because it would add complexity to the model since an additional index would be needed that contemplates the temporal component. However, the application presented below does take into account this temporal component. Once we provide as input a set of orders, the output is all the daily routes needed to fulfil orders completely. However, this application, for each day, is also governed to maximize the load, taking care of the orders requested for that day and without taking into account the urgency of the remaining orders. 4.4 Examples To illustrate the model that we defined, we consider two small instances of data. 4.4.1 A fictitious example In the first example, we chose to work with a fleet of two trucks: a small truck with two hoppers that can each hold up to 4 and 3 of feeds and a larger truck with three hoppers, whose capacities are 5, 4 and 4 of feed. The first truck cannot carry more than 6, while the other can carry up to 20 at a time. To reduce the number of iterations, let us suppose that each truck can only follow one route per day. Moreover, there are six stockbreeders located in Orense, Begonte, Friol, Sarria, A Estrada and Forcarei. We assume that the larger truck is unable to deliver food to the stockbreeders located in Begonte and Sarria. In addition, to simplify the problem and to make the resolution faster, we suppose that the orders are always for only one type of feed (this assumption is well justified since it corresponds to what the stockbreeders usually ask for) and that the amount of feed in each order is close to the capacity of one of the hoppers. The orders are as shown in Table 1 , and the distances between the different barns and the times required to go from one to another are as listed in Table 2. Table 1. Orders of different stockbreeders (in tonnes of feed) Orders Feed 1 Feed 2 Feed 3 Feed 4 Days left to deliver Orense 4.5 1 Begonte 4 1 Friol 3 1 Sarria 4 3 A Estrada 3.5 1 Forcarei 2.5 2 Orders Feed 1 Feed 2 Feed 3 Feed 4 Days left to deliver Orense 4.5 1 Begonte 4 1 Friol 3 1 Sarria 4 3 A Estrada 3.5 1 Forcarei 2.5 2 Open in new tab Table 1. Orders of different stockbreeders (in tonnes of feed) Orders Feed 1 Feed 2 Feed 3 Feed 4 Days left to deliver Orense 4.5 1 Begonte 4 1 Friol 3 1 Sarria 4 3 A Estrada 3.5 1 Forcarei 2.5 2 Orders Feed 1 Feed 2 Feed 3 Feed 4 Days left to deliver Orense 4.5 1 Begonte 4 1 Friol 3 1 Sarria 4 3 A Estrada 3.5 1 Forcarei 2.5 2 Open in new tab Table 2. Distances and traveling times between pairs of stockbreeders km/Min Coop Orense Begonte Friol Sarria A Estrada Forcarei Coop. 0 50/47 61/60 44/41 46/46 85/69 70/67 Orense 0 111/104 93/86 78/77 96/62 81/60 Begonte 0 19/20 56/44 116/97 156/105 Friol 0 57/55 103/90 92/89 Sarria 0 144/115 129/113 A Estrada 0 20/17 Forcarei 0 km/Min Coop Orense Begonte Friol Sarria A Estrada Forcarei Coop. 0 50/47 61/60 44/41 46/46 85/69 70/67 Orense 0 111/104 93/86 78/77 96/62 81/60 Begonte 0 19/20 56/44 116/97 156/105 Friol 0 57/55 103/90 92/89 Sarria 0 144/115 129/113 A Estrada 0 20/17 Forcarei 0 Open in new tab Table 2. Distances and traveling times between pairs of stockbreeders km/Min Coop Orense Begonte Friol Sarria A Estrada Forcarei Coop. 0 50/47 61/60 44/41 46/46 85/69 70/67 Orense 0 111/104 93/86 78/77 96/62 81/60 Begonte 0 19/20 56/44 116/97 156/105 Friol 0 57/55 103/90 92/89 Sarria 0 144/115 129/113 A Estrada 0 20/17 Forcarei 0 km/Min Coop Orense Begonte Friol Sarria A Estrada Forcarei Coop. 0 50/47 61/60 44/41 46/46 85/69 70/67 Orense 0 111/104 93/86 78/77 96/62 81/60 Begonte 0 19/20 56/44 116/97 156/105 Friol 0 57/55 103/90 92/89 Sarria 0 144/115 129/113 A Estrada 0 20/17 Forcarei 0 Open in new tab We define the values of the other data as follows. As stated before, all the costs given here are arbitrary and may not actually be relevant. Maximum daily distance covered by a truck: 600 km. Cost of an unloading: €2 per delivery. Fixed cost of food transportation: €1 per tonne carried. Variable cost of food transportation: – For 0 to 100 km covered, €0.5 per km per tonne carried. – For 100 to 200 km covered, 0.75 euros per km per tonne carried. – For 200 to 300 km covered, €1 per km per tonne carried. The model was implemented in the AMPL language (cf. Fourer et al., 2003). We solved the problem using a free version of the solver KNITRO (cf. Byrd et al., 2006), which uses a version of the interior point algorithm for non-linear problems, obtained in the NEOS optimization server. When dealing with integer variables, the solver uses a branch and bound algorithm that can only give a locally optimal solution. For maximization, we configured the solvers to try multiple starting points, 50 to be exact, to have the result be as accurate as possible. KNITRO manages to find a locally optimal solution for 17 in a total processing time of 1180 s (approximately 20 min). Here, we are not able to satisfy all the non-urgent orders. Indeed, we can bring the entire order of the stockbreeder located in Friol, but only 2 out of 4 to the stockbreeder in Sarria and nothing to the stockbreeder in Forcarei. Note that it is feasible to serve a client on different days. This is not a problem for the customer and can contribute to the company’s goal of maximizing the daily load. We can now try to find a satisfactory solution to the whole problem, again with 50 different starting points. The locally optimal solution that the solvers finally give after 9278 s (approximately 2.5 h) is as follows. The larger truck travels from Taboada to Orense, then A Estrada and finally Friol and carries – 4.5 of feed 3 to the stockbreeder in Orense in its first hopper; – 3.5 of feed 4 to the stockbreeder in A Estrada in its second hopper; – 3 of feed 2 to the stockbreeder in Friol in its third hopper. The small truck goes from Taboada to Sarria and then Begonte and carries – 4 of feed 1 to the stockbreeder in Begonte in its first hopper; – 2 of feed 4 to the stockbreeder in Sarria in its second hopper. In this configuration, the total cost of the deliveries is €3966.5. 4.4.2 An example with real data In this second example, we consider real data provided by the cooperative AIRA. We take two similar trucks with five hoppers that can contain up to 3, 3.7, 3.8, 3.7 and 3 of feed. The trucks cannot carry more than 11.6 (thus, one truck is not able to fulfil all of the orders). Let us suppose that the trucks can follow one route per day. Moreover, there are five stockbreeders, all of which have urgent orders of only one type of feed. Each truck can deliver food to all the stockbreeders. In addition, to simplify the problem and make resolution faster (and also because all the orders are urgent), we suppose that the objective is only to minimize the total distance travelled. Note that the final number of hoppers used is larger than the number of served stockbreeders. The solution shows how we have to share each order among different hoppers. In the previous example, each fulfilled order was allocated to only one hopper because of the features of the problem. The orders are as shown in Table 3 , and the distances between the different barns are as given in Table 4. Table 3. Orders of the different stockbreeders (in tonnes of feed) Orders Feed 1 Days left to deliver Stockbreeder 1 3.300 1 Stockbreeder 2 2.951 1 Stockbreeder 3 3.003 1 Stockbreeder 4 3.016 1 Stockbreeder 5 2.496 1 Orders Feed 1 Days left to deliver Stockbreeder 1 3.300 1 Stockbreeder 2 2.951 1 Stockbreeder 3 3.003 1 Stockbreeder 4 3.016 1 Stockbreeder 5 2.496 1 Open in new tab Table 3. Orders of the different stockbreeders (in tonnes of feed) Orders Feed 1 Days left to deliver Stockbreeder 1 3.300 1 Stockbreeder 2 2.951 1 Stockbreeder 3 3.003 1 Stockbreeder 4 3.016 1 Stockbreeder 5 2.496 1 Orders Feed 1 Days left to deliver Stockbreeder 1 3.300 1 Stockbreeder 2 2.951 1 Stockbreeder 3 3.003 1 Stockbreeder 4 3.016 1 Stockbreeder 5 2.496 1 Open in new tab Table 4. Travelling distances between pairs of stockbreeders (in km) Min Coop Stockb. 1 Stockb. 2 Stockb. 3 Stockb. 4 Stockb. 5 Coop. 0 28 69 64 27 17 Stockb. 1 0 67 62 20 20 Stockb. 2 0 7 74 58 Stockb. 3 0 69 53 Stockb. 4 0 25 Stockb. 5 0 Min Coop Stockb. 1 Stockb. 2 Stockb. 3 Stockb. 4 Stockb. 5 Coop. 0 28 69 64 27 17 Stockb. 1 0 67 62 20 20 Stockb. 2 0 7 74 58 Stockb. 3 0 69 53 Stockb. 4 0 25 Stockb. 5 0 Open in new tab Table 4. Travelling distances between pairs of stockbreeders (in km) Min Coop Stockb. 1 Stockb. 2 Stockb. 3 Stockb. 4 Stockb. 5 Coop. 0 28 69 64 27 17 Stockb. 1 0 67 62 20 20 Stockb. 2 0 7 74 58 Stockb. 3 0 69 53 Stockb. 4 0 25 Stockb. 5 0 Min Coop Stockb. 1 Stockb. 2 Stockb. 3 Stockb. 4 Stockb. 5 Coop. 0 28 69 64 27 17 Stockb. 1 0 67 62 20 20 Stockb. 2 0 7 74 58 Stockb. 3 0 69 53 Stockb. 4 0 25 Stockb. 5 0 Open in new tab The locally optimal solution that the solvers provide after 109.92 s is as follows. Truck 1 goes from Taboada to stockbreeder 5, then to stockbreeder 3 and then to stockbreeder 2 and carries – 1.47 of feed (49.8|$\%$| of the demand) to stockbreeder 2 in hopper 1, – 2.50 of feed (100|$\%$| of the demand) to stockbreeder 5 in hopper 2, – 1.55 of feed (51.6|$\%$| of the demand) to stockbreeder 3 in hopper 3, – 1.45 of feed (48.4|$\%$| of the demand) to stockbreeder 3 in hopper 4 and – 1.47 of feed (50.2|$\%$| of the demand) to stockbreeder 2 in hopper 5. Truck 2 goes from Taboada to stockbreeder 4 and then to stockbreeder 1 and carries – 1.51 of feed (50|$\%$| of the demand) to stockbreeder 4 in hopper 1, – 3.30 of feed (100|$\%$| of the demand) to stockbreeder 1 in hopper 3 and – 1.51 of feed (50|$\%$| of the demand) to stockbreeder 4 in hopper 5. In this configuration, the total distance of the deliveries is 221 km. 5. Routing implementation through heuristic algorithms Because our problem is hard to solve exactly, we design and apply a hybrid heuristic approach that combines an insertion algorithm to obtain an initial solution and then uses a simulated annealing metaheuristic to improve it, because these kinds of strategies have been shown to be effective tools for solving hard combinatorial optimization problems. 5.1 A heuristic algorithm for obtaining an initial solution This algorithm is created with the aim of calculating an initial solution to the problem in a short run time (seconds). While this solution is not optimal, the method requires a quality result from which to iterate future optimizations. We start with the implementation of the data structure of the problem, with all the attributes defined in the statement provided by the cooperative. In this implementation, the algorithm progressively develops truck routes that serve customers until all the orders have been planned. The pseudocode for this algorithm is explained below. Figure 3 illustrates the algorithm. The programming of each journey begins with the choice of a customer to be served (client seed) and an initial truck (lorry seed). If there are no customers to serve, the algorithm terminates. If no truck can serve a customer, it is understood that the daily work is planned for all trucks and the program is restarted (with empty vehicles) the next day to serve the pending customers. If programming has been completed for 365 days and there are still customers who have not been served, the algorithm terminates. Build the path from the central repository, visit the customer seed and return to the company. Load truck hoppers in the most efficient way. If no additional customer can be inserted in this step, the truck has completed its daily schedule; return to step 1 to define more routes for other trucks. Search all the customers who have not been served and their chances of insertion into each subpath of the route. A client and a position are chosen such that the lowest cost in the insertion is ensured. The selected customer, along with the position, is inserted into the journey. The hoppers are loaded in all subpaths that precede this position. Check that the restrictions of the problem are not violated (if they are violated, return to step 3 and discard this client). In the case that an order does not fit completely in the hoppers of the assigned vehicle, the request is divided into two: one portion is loaded into the free space of the truck and is considered to be fulfilled, and the other is not. Return to step 2 to continue inserting customers. If in step 2 no other client can be inserted, we consider the route to be finalized. Return to step 1 to define more routes with additional trucks. Fig. 3. Open in new tabDownload slide Diagram of the insertion algorithm used. Fig. 3. Open in new tabDownload slide Diagram of the insertion algorithm used. The user can select the strategy of customer seed choice from the following options: farthest from the factory, customer with the highest workload or a random customer. The user can also select the strategy for choosing the truck to use in each case from the following options: the truck with the lowest mileage, the highest capacity truck or a random truck. In this way, we can implement nine types of initial solution according to the selection of the seed customer and the seed truck: 1i Customer farthest from the factory and the truck with the lowest mileage, 2i Customer farthest from the factory and the highest capacity truck, 3i Customer farthest from the factory and a random truck. 4i Customer with the highest workload and the truck with the lowest mileage, 5i Customer with the highest workload and the highest capacity truck, 6i Customer with the highest workload and a random truck, 7i A random customer and the truck with the lowest mileage, 8i A random customer and the highest capacity truck and 9i A random customer and a random truck. 5.2 The simulated annealing algorithm The algorithm acts on the initial solution, looking for possible improvements that are managed and optimized in subsequent iterations. The structure of this methodology is as follows. While the maximum number of iterations or time has not been exceeded, the algorithm will try to find a neighbouring solution to the current one using random methods. Execution ends if the maximum number of iterations or time allowed is reached or if it is not possible to find any neighbouring solution. In that case, the lowest-cost processing result is returned. At each iteration, the search process moves from the current solution to a neighbouring solution. According to the classical simulated annealing approach, let |$T$| be a parameter that measures the tendency to accept the current candidate. Here |$T$| is called the temperature. Among the immediate neighbours of the current solution, select one randomly to become a potential future solution. If the candidate under consideration is better than the current solution, it is accepted. If the candidate under consideration is worse, the chance of acceptance will depend on how much worse it is and on the temperature value, which will evolve over the execution of the algorithm (it will decrease, thus decreasing the associated probability). The value of |$T$| is updated, if applicable, and we return to step 1. For more details about simulated annealing algorithms, see Nunes De Castro (2006). To implement an efficient software, seven ways to obtain neighbouring solutions are designed and programmed and are randomly interspersed. Only in this way will we be able to explore a sufficiently large part of the region of potential solutions, where the global optimum can be located, if the execution time permits. Next, the implemented procedures are described. Figure 4 illustrates the algorithm. 1n Obtain a new neighbour by changing the orders of a randomly selected day, truck, customer and journey to a different randomly selected truck, journey and day. 2n Obtain a new neighbour by moving the orders of one randomly selected day, truck, journey and client to a different journey of the same truck and day. 3n Obtain a new neighbour by moving the orders of a day, truck, customer and journey to a random truck and journey of the same day. 4n Obtain a new neighbour by exchanging the orders of a day, truck, journey and random customer with those of another random customer and journey of the same truck and day. 5n Obtain a new neighbour by exchanging the orders of a randomly selected day, truck, customer and journey with those of another randomly selected different day, truck, customer and journey. 6n Obtain a new neighbour by randomly choosing 1 day, a truck scheduled for that day, a path of that truck and two clients served in that path. We exchange the order of these customers for that truck, route and day. 7n Obtain a new neighbour by randomly choosing 1 day, two trucks scheduled for that day, a path for each truck and one client served in each path. We exchange the farmer served in a path from one truck to the other on the other truck. Fig. 4. Open in new tabDownload slide Diagram of the simulated annealing algorithm used. Fig. 4. Open in new tabDownload slide Diagram of the simulated annealing algorithm used. Simulated annealing performance usually depends on the calibration of its parameters, namely the number of iterations at the same temperature and the law for temperature decrease at each macro iteration. In our experiments, during the execution of the algorithm, the temperature was updated 10 times, although the user could define the number of total resources (either iterations or time) that the entire algorithm used. For example, if we executed the algorithm specifying a maximum of 1000 iterations or 20 s, the temperature would be updated upon reaching 100 iterations or at 2 s, then upon reaching 200 iterations or at 4 s and so on. The temperature starts from a value proportional to the cost obtained by the initial solution (the heuristic insertion algorithm), specifically this cost multiplied by 0.9. Each time the temperature is updated, it is proportionally reduced by multiplying the current temperature by 0.6. 6. Some experimental results This section describes the results obtained from real and simulated examples and presents evidence to verify that the requirements of our application are satisfied. In most of these examples, more complex situations are contemplated than in the examples seen so far, in the sense that routes are contemplated on several days and trucks have to travel more than one route per day. To perform the computations, a personal computer was used with a processor Intel Core i5 2.5 GHz and 6.00 GB of RAM. 6.1 A comparison between the optimal solution and that calculated with the heuristic algorithms We return to the example with real data from Section 4.4.2 and we solve it using our algorithms. In Table 5, we note the distance travelled that originates the solution obtained by the insertion algorithm in the different possible cases of choice of seed client and truck to be used. Table 5. Results of the heuristic insertion algorithm (expressed in km) for the example with real data Truck/client Farthest More urgent Random Less distance 247 221 221 More capacity 221 319 247 Random 247 319 247 Truck/client Farthest More urgent Random Less distance 247 221 221 More capacity 221 319 247 Random 247 319 247 Open in new tab Table 5. Results of the heuristic insertion algorithm (expressed in km) for the example with real data Truck/client Farthest More urgent Random Less distance 247 221 221 More capacity 221 319 247 Random 247 319 247 Truck/client Farthest More urgent Random Less distance 247 221 221 More capacity 221 319 247 Random 247 319 247 Open in new tab With respect to the algorithm of simulated annealing, in all cases a distance of 221 km is obtained. As we saw in Section 4.4.2, the optimal solution provided by KNITRO is 221 km, which coincides with that obtained through the simulated annealing algorithm we have designed. It should also be mentioned that although the distance obtained is the same, the solution is not identical. In particular, with KNITRO we obtained that the orders of some customers were transported in different hoppers, while the solution provided by our algorithm, in this case, is such that the order of each customer is transported in a single hopper. In addition, here we allow only one truck to be used, which travels two routes. The initial solution obtained by the insertion algorithm also obtains the value of 221 km in some cases, although in others it provides other local optima. In particular, we see in Table 5 that the best results of the insertion algorithm are obtained, in this example, when the nearest truck is used and the customer is selected at random or when the farthest customer is selected and the truck with greatest capacity is taken. 6.2 Tests with real data on results in terms of costs In collaboration with the agricultural cooperative, we were able to obtain information about our problem, and so to generate examples of real instances. In this section, we work with one of the real examples to illustrate the operating characteristics of the software. In this case, we pay special attention to the costs associated with the routes designed. The example involves the company located in the municipality of Taboada with a fleet of two trucks, each with five hoppers. Its customers consist of 21 farmers scattered across Galicia. Customers place a total of 22 urgent requests for a single feed. Files are provided with the distance matrix, load characteristics and hoppers of the trucks and data from the orders. The information could correspond to the planning of a day of work of the company. To enter all the data of the problem into the software, it is necessary to generate the required entities (trucks, locations, orders etc.) and define their characteristics using the graphical interface. After this process is completed, we are ready to run the programmed algorithms with the desired characteristics. First, the insertion heuristic methodology is executed to obtain initial solutions. This algorithm operates under two different parameters: the form of the seed customer choice and the initial truck. Table 6 lists the values of the objective function in Euros resulting from this execution, along with each of the parameters of the function. Table 6. Real results of the heuristic insertion algorithm (expressed in Euros) Truck/client Farthest More urgent Random Less distance 1008 1008 1029 More capacity 1007 1007 1054 Random 1005 1010 1040 Truck/client Farthest More urgent Random Less distance 1008 1008 1029 More capacity 1007 1007 1054 Random 1005 1010 1040 Open in new tab Table 6. Real results of the heuristic insertion algorithm (expressed in Euros) Truck/client Farthest More urgent Random Less distance 1008 1008 1029 More capacity 1007 1007 1054 Random 1005 1010 1040 Truck/client Farthest More urgent Random Less distance 1008 1008 1029 More capacity 1007 1007 1054 Random 1005 1010 1040 Open in new tab The different initial solutions offer a minimum score of €1005 and an average result of €1019. Subsequently, the algorithm is executed on this data to implement the simulated annealing metaheuristic methodology to search for better solutions. For each situation, we have executed a total of 50000 iterations or have employed 300-s runs, with the limiting factor being the time in most cases. Table 7 shows the solutions improved by this process for each of the above results. Table 7. Real results of the simulated annealing algorithm metaheuristic (expressed in Euros) Truck/client Farthest More urgent Random Less distance 788 775 838 More capacity 921 992 848 Random 917 992 842 Truck/client Farthest More urgent Random Less distance 788 775 838 More capacity 921 992 848 Random 917 992 842 Open in new tab Table 7. Real results of the simulated annealing algorithm metaheuristic (expressed in Euros) Truck/client Farthest More urgent Random Less distance 788 775 838 More capacity 921 992 848 Random 917 992 842 Truck/client Farthest More urgent Random Less distance 788 775 838 More capacity 921 992 848 Random 917 992 842 Open in new tab After the optimization process is completed, the solutions offer a minimal result of €775 and an average result of 879 euros. We note that the first two results in the first row of Table 7 may seem a little surprising. We believe that a possible explanation is that the initial solutions from which they depart, although in terms of objective function are similar to all others, are not so much in terms of variables and therefore of position in the feasible region and, hence, they can lead to a better result after execution of the simulated annealing algorithm. In any case, we obtain an average saving of 13.8|$\%$|⁠, which can reach 23|$\%$| in the best cases. These are very interesting economic values for the case of a transport logistics activity. In view of the results, we can highlight the homogeneity of the initial solutions, among which there is very little variability, at least in this example. This strengthens the idea that with any selection method used to obtain the seed customers or initial truck the solutions that are produced are of similar quality after the instant execution. The operation of the optimization algorithm for a moderate length of time (5 min) generates improvements in all these situations, and many of these improvements involve substantial amounts of money. This justifies, in principle, the use of this type of methodology in resolving VRPs. Fig. 5. Open in new tabDownload slide Improvement between initial and optimized solution on problems of different sizes of data set. Fig. 5. Open in new tabDownload slide Improvement between initial and optimized solution on problems of different sizes of data set. 6.3 Tests on the efficiency of the application With the numerical results of different simulated examples, we have constructed two figures to compare the degree of improvement in the results provided by the simulated annealing algorithm with the data set size and the runtime of the algorithm. From analysis of the figures, we can obtain information about the efficiency of the application. Figure 5 shows how the system behaves as it faces problems of different sizes. The vertical axis indicates the percentage of improvement between the initial solution and optimized solution, while the horizontal axis shows the size of the data set. The data set includes the distance matrix, load characteristics and hoppers of the trucks and the orders. The optimization runs for a constant time of 30 s in each case. Fig. 6. Open in new tabDownload slide Improvement between initial and optimized solution depending on the time needed to generate them. Fig. 6. Open in new tabDownload slide Improvement between initial and optimized solution depending on the time needed to generate them. Figure 5 shows a function similar to an increasing linear function. In fact, the value of the linear correlation coefficient is 0.93, which gives rise to a coefficient of determination of 0.87. The proximity of the graph to the origin of the coordinates indicates that in very small problems, the optimization is hardly effective. This supports the hypothesis that the results of the heuristic algorithm of insertion are of high quality. In uncomplicated cases, the initial solution provided by this algorithm cannot be improved significantly, so it is already very close to the optimum. As the problem grows in size, it becomes impossible to obtain optimal solutions with heuristic methods using instant execution. It is therefore effective to use the metaheuristic optimization algorithms, which obtain better solutions during the allowed execution time. As we can see in the figure, the simulated annealing algorithm performs efficiently on problems of relatively large sizes. Figure 6 shows, in a more specific way, the concept of efficiency. It measures the quality of a result as a function of the time necessary to generate it. It indicates how the results of the optimization process evolve as the runtime changes. The vertical axis corresponds to the percentage of improvement between the initial solution and optimized solution, and the horizontal axis corresponds to the time spent obtaining the result. In all these tests, the same problem of moderate size was used. The increasing character of the function tells us that longer runtimes produce higher-quality optimizations. If we deepen the analysis, it highlights the proximity of the graph to the origin of the coordinates, indicating that with very small execution times, there are no significant improvements. In those cases, the optimization algorithm does not have enough processing time to find a local minimum of the objective function to improve the initial solution. As the execution time increases, usually new local minima that improve the initial solution emerge and, progressively, the identification of the best solutions is optimized as the number of results increases. Then, the rate of improvement increases and the graph reaches its maximum slope. It is in this position that the temporary resource is exploited and the appearance of higher-quality results is more frequent. When the runtime becomes large, the results improve, but at decreasing rates. We believe that in these cases, the solution obtained is close to the global optimum of the objective function. The refinement process of the results can continue, but the quality of the solutions will not improve significantly. In view of the obtained results, tests and subsequent analysis, the initial requirements for which this application was designed are validated and accepted. 6.4 Planning for more than one day and post-optimality analysis In Sections 6.4 and 6.5, we undertake a more detailed examination of the results corresponding to planning the work of 2 days or more, by assuming that each truck can travel more than one route each day. In this section, we used real data in which we have 36 orders, 19 of them being urgent. The total amount of feed requested is 205, with the most distant order being 65 km away. We ran the simulated annealing algorithm once with each of the possible initial solutions. We first considered two trucks with a total capacity of 30.6. In Fig. 7 (dotted line), we show the resulting distances travelled. The horizontal axis contains nine different types of initial solution (i.e. taking into account the three possible strategies of customer seed choice, and for each of them, the three strategies for choosing the truck). The vertical axis shows the distance corresponding to each solution. In this case, we needed 2 days to service all orders. We also performed a sensitivity analysis or post-optimality analysis by varying the dimensions of the fleet of trucks available and we also collected the results in Fig. 7. Following this analysis, the solid line refers to the results when we doubled the number of trucks; that is, we now have a capacity of 61.2, and the dashed line was constructed after considering three larger trucks with a total capacity of 54. In the two latter cases, 1 day is sufficient to cover all orders. It is obvious that, overall, the best results were obtained with three large trucks. Within each series, the best result is obtained when the seed client is randomly selected in the initial solution. Fig. 7. Open in new tabDownload slide Comparative of the simultaneous planning of 36 orders corresponding to 2 real working days using different sizes of the truck fleet. Fig. 7. Open in new tabDownload slide Comparative of the simultaneous planning of 36 orders corresponding to 2 real working days using different sizes of the truck fleet. 6.5 Variability of results when we are planning for more than one day In this experiment, we consider a simulated data set corresponding to a company with three trucks with a total capacity of 40. We have 65 clients and 59 orders, with a total requested feed of 561. Only 15 orders are urgent and we have an average delivery time of 2.15 days. The maximum distance from the factory to a customer is 62 km. We run the simulated annealing algorithm 10 times, for 15 min at a time, with each of the nine different methods of obtaining the initial solution. In all cases, the results give rise to 3 days to attend all orders, always respecting the delivery time. On the first 2 days, all the trucks travel almost all permitted kilometres, making use of almost the entire load capacity, and every day at least one truck travels more than one route. With the results obtained, in terms of distances expressed in kilometres, we constructed Fig. 8, which contains nine box plots, one for each procedure of obtaining the initial solution (indicated on the horizontal axis) and taking as a variable the total distance travelled in kilometres (indicated on the vertical axis). The interest of this graph lies in comparing once again the quality of the different solutions that our application can provide and, in this case, also the variability of each of them. As we see, in this example, there are six procedures whose best solution is very close. Among them, the type of solution that provides the best solution is also the one that presents results with more variability. Among the different procedures, we see that there are differences in the variability and the symmetry of the result sets. If we look at the medians (thick line), the best results (and also one of the worst ones) correspond to cases with null variability. In view of the results, as we are dealing with executions of 15 min, it seems advisable, before each planning, to make several runs considering the different procedures and select the most favourable. In addition, as seen in the previous section, analysis of the effect of the use of different types of trucks can be useful. Fig. 8. Open in new tabDownload slide Simultaneous planning of 3 simulated working days (variability of the obtained results when running the algorithm 10 times). Fig. 8. Open in new tabDownload slide Simultaneous planning of 3 simulated working days (variability of the obtained results when running the algorithm 10 times). 6.6 Effectiveness of different types of neighbourhood In the simulated annealing algorithm designed, we have seen that seven different types of neighbourhood structures were used randomly. The probabilities with which each of these types of neighbourhood are chosen are different in order to achieve a better functioning and as advised different tests performed. Different experiments have been carried out to explore the efficacy of each neighbourhood structure separately and we have come to the conclusion that such efficacy depends largely on the particular problem to be solved, so it seems a good strategy to take all of them into account for getting a more efficient algorithm. Figure 9 shows again the example considered in Section 6.4, making use of the three large trucks. Now, we consider separately the seven different neighbourhood structures (1n–7n) presented in Section 5.2, and in addition, two other neighbourhood structures combinations: 8n A combination with equal probability of the seven neighbourhood structures 1n–7n. 9n Our original algorithm that uses the seven neighbourhood structures 1n–7n with different probabilities. Fig. 9. Open in new tabDownload slide Simultaneous planning of 36 orders corresponding to 2 real working days, using three large trucks (variability of the obtained results when running the algorithm 30 times). Fig. 9. Open in new tabDownload slide Simultaneous planning of 36 orders corresponding to 2 real working days, using three large trucks (variability of the obtained results when running the algorithm 30 times). In all cases, the same procedure is used to obtain the initial solution. In the vertical axis, we represent the distances travelled as a result of our planning. We run each of the algorithms 30 times for 15 minutes of time. In this case, neighbourhood structures 1n, 2n and 7n are the ones that offer better results on average, although in some cases with more dispersion than the algorithms that use all of them. In addition, with other data, it is not always these three structures that offer the best results. In any case, in the last version of our interface we have incorporated the possibility of choosing among the different neighbourhood structures when executing the algorithm. This can be useful for an expert user, to carry out research, or if we think about a system export to other companies. 6.7 A longer-term planning In this last section of results, we show the capacity of the application to plan the transport of the feed factory for a longer period of time, through simulated data. As for the orders to be served, they constitute a number of 720 and represent a total of 3229. The furthest order from the factory is 57 km away. Each of the orders is between 4 and 6, and in terms of urgency, it is very relaxed and everyone can be served any day of work for the 2 weeks following the time of planning. With respect to the truck fleet, there are two trucks, with a total capacity of 60, five hoppers each truck and with a restriction of at most 600 km per truck and day. After 5400 s running the application, we obtain a planning that allows us to meet the orders in 7 working days. The main characteristics of the planning are shown in Table 8. Table 8. Kilometres, number of routes made by trucks and load in the planning of the distribution of 720 orders in 7 days Day 1 2 3 4 5 6 7 Km truck 1 595 597 593 598 600 592 585 Number of routes truck 1 19 10 11 11 10 8 6 Km truck 2 594 600 595 597 599 575 298 Number of routes truck 2 15 13 9 9 11 8 5 Average load 23|$\%$| 31|$\%$| 32|$\%$| 29|$\%$| 32|$\%$| 30|$\%$| 34|$\%$| Day 1 2 3 4 5 6 7 Km truck 1 595 597 593 598 600 592 585 Number of routes truck 1 19 10 11 11 10 8 6 Km truck 2 594 600 595 597 599 575 298 Number of routes truck 2 15 13 9 9 11 8 5 Average load 23|$\%$| 31|$\%$| 32|$\%$| 29|$\%$| 32|$\%$| 30|$\%$| 34|$\%$| Open in new tab Table 8. Kilometres, number of routes made by trucks and load in the planning of the distribution of 720 orders in 7 days Day 1 2 3 4 5 6 7 Km truck 1 595 597 593 598 600 592 585 Number of routes truck 1 19 10 11 11 10 8 6 Km truck 2 594 600 595 597 599 575 298 Number of routes truck 2 15 13 9 9 11 8 5 Average load 23|$\%$| 31|$\%$| 32|$\%$| 29|$\%$| 32|$\%$| 30|$\%$| 34|$\%$| Day 1 2 3 4 5 6 7 Km truck 1 595 597 593 598 600 592 585 Number of routes truck 1 19 10 11 11 10 8 6 Km truck 2 594 600 595 597 599 575 298 Number of routes truck 2 15 13 9 9 11 8 5 Average load 23|$\%$| 31|$\%$| 32|$\%$| 29|$\%$| 32|$\%$| 30|$\%$| 34|$\%$| Open in new tab As we can see, both trucks drive very close to the maximum number of kilometres allowed every day, except one of them on the last day. As for the number of routes, it shows a decreasing trend with the passage of days. In this way, it is possible to advance the deliveries of the orders. Regarding the average load, it varies between 23|$\%$| and 34|$\%$|⁠, which is due to the restriction that different types of feed or orders of different customers can not be mixed in the same hopper. 7. Managerial implications Among the many problems raised by the transfer of mathematics to industry, we have investigated the efficient transport of food carried out by many companies such as agricultural cooperatives that need to supply animal food to farms. In Galicia, the Spanish region where the cooperative that has inspired this work is located, we can find more than 15 companies that manufacture feed and in Spain the number is more than 500. The contribution of the techniques of operations research in these and other logistical problems is becoming increasingly fruitful and this is a fact that is reflected in the scientific literature, both theoretical and applied. But even now, some small and medium businesses make use of a ‘human-based’ methodology in transport management, which requires much time to be devoted to such management and sometimes with inefficient results. On other occasions, commercial tools that are expensive are used, and as they have not been designed for a specific company, they are laborious to adapt and they can be a little transparent for the end users. We believe that the tool here proposed facilitates the task of managing truck routes by the logistics technician of the cooperative. Another positive aspect of our decision support tool, which we feel is appropriate to point out, is the possibility of extending it to similar companies. A situation very close to the one we have studied is the transport of cereals. Another example, which in principle would not require major modification, is the dairy factories that need to transport milk from the farms to the factory. Another possible extension would be the transport of biofuel. To cite one more example, also agricultural and related to livestock, this would be the transport of animals from farms to a slaughterhouse. In the latter case, it would be necessary to adapt the set of technical constraints that model the unloading process and probably the objective cost function. In any case, when applied to any extension, it is necessary to adjust parameters, both of the model itself and of those related to execution, such as number of iterations or advisable time. 8. Conclusions and final remarks We have proposed a tool designed according to the specific requirements of a company. We can conclude that our first experimental results show that the designed algorithms can be applied effectively with reasonable computational effort. In addition, we have created an interface that allows one to use the algorithms in a friendly way that can make the actual implementation of our methodology possible in the cooperative. It is worth noting the ease with which post-optimality analysis can be performed with our tool, such as the related with the selection of trucks to be used within the available fleet, which may affect the result of the designed route and the associated costs. Furthermore, the comparison between our results and those provided by the former methodology shows a reduction both in working time and economic cost. Thus, the advantage that the use of this tool can provide is a reduction in distribution costs because the distances obtained through an optimization procedure will be lower than those obtained manually without using the program. In the experiments we have included in the work, regarding schedules beyond 1 day, we saw how the application is able to schedule deliveries, maximizing the daily load and respecting the requests of the partners. It is important to foster a culture in the company that allows reliable forecasts, carrying out a rigourous control of orders and deliveries, remnants at the time of ordering and that the partners make the orders with reasonable notice and all relevant incidents are communicated to the company, to be able to avoid saturations of the system when fulfilling the orders. Moreover, this tool can easily connect with others, either aimed at automating the various stages of transport and other more administrative tasks such as report generation. On the other hand, there are other metaheuristics defined for problems of vehicle routes so that future research could include comparing the result of its application with that obtained in this work. In particular, for an application in a longer-range planning, it may be of interest to deepen the behaviour of the different neighbourhood structures considered in the simulated annealing algorithm, and rather than to make a random selection in the same, follow some type of strategy, as well as perform path relinking to intensify the search after finding local optima. Finally, we can say that along with the distribution problem studied here, there are other aspects associated with possible inefficiencies in the system that might be worth studying (cf. Amiama-Ares et al., 2014), such as the problem of the production of feed, which can be addressed in the future using the tools of operations research and computer science. Acknowledgements We appreciate the comments of M. Álvarez, C. Amiama, L. Carpente, M. Mosquera, M. C. Perrot, an associate editor and two anonymous referees on this issue. Funding Ministerio de Economía y Competitividad and European Regional Development Fund (ERDF) (MTM2014-53395-C3-2-P, MTM2017-87197-C3-3-P). References Ambrosino , D. & Sciomachen , A. ( 2007 ) A food distribution network problem: a case study . IMA J. Manage. Math. , 18 , 33 – 53 . Google Scholar Crossref Search ADS WorldCat Amiama-Ares , C. , Loza-García , J. R. & Bueno-Lema , J. ( 2014 ) Spatial decision support system for the route management of feed delivery . 18th International Congress on Project Management and Engineering . Alcañiz, Teruel, Spain : University of Zaragoza and the Spanish Association of Project Management and Engineering . Google Preview WorldCat COPAC Bent , R. & van Hentenryck , P. ( 2004 ) A two-stage hybrid local search for the vehicle routing problem with time windows . Transport. Sci. , 38 , 515 – 530 . Google Scholar Crossref Search ADS WorldCat Bouzaiene-Ayari , B. , Dror , M. & Laporte , G. ( 1993 ) Vehicle routing with stochastic demands and split deliveries . F. C. D. S. , 18 , 63 – 69 . WorldCat Brito , J. , Expósito , A. & Moreno , J. ( 2016 ) Variable neighbourhood search for close-open vehicle routing problem with time windows . IMA J. Manage. Math. , 27 , 25 – 38 . Google Scholar Crossref Search ADS WorldCat Byrd , R. H. , Nocedal , J. & Waltz , R. A. ( 2006 ) KNITRO: an integrated packages for nonlinear optimization . Large-Scale Nonlinear Optimization ( G. di Pillo & M. Roma eds ). Boston : Springer . Google Preview WorldCat COPAC Carpente , L. , Casas-Méndez , B. , Jácome , C. & Puerto , J. ( 2010 ) A model and two heuristic approaches for a forage harvester planning problem: a case study . TOP , 18 , 122 – 139 . Google Scholar Crossref Search ADS WorldCat Chao , I. M. ( 2002 ) A tabu search method for the truck and trailer routing problem . Comput. Oper. Res. , 29 , 33 – 51 . Google Scholar Crossref Search ADS WorldCat Chen , S. , Golden , B. & Wasil , E. ( 2007 ) The split delivery vehicle routing problem: applications, algorithms, test problems, and computational results . Networks , 49 , 318 – 329 . Google Scholar Crossref Search ADS WorldCat Clarke , G. & Wright , W. ( 1964 ) Scheduling of vehicles from a central depot to a number of delivery points . Oper. Res. , 12 , 568 – 581 . Google Scholar Crossref Search ADS WorldCat Dantzig , G. B. & Ramser , J.H . ( 1959 ) The truck dispatching problem . Manag. Sci. , 6 , 80 – 91 . Google Scholar Crossref Search ADS WorldCat Derigs , U. , Gottlieb , J. , Kalkoff , J. , Piesche , M. , Rothlauf , F. & Vogel , U. ( 2011 ) Vehicle routing with compartments: applications, modelling and heuristics . OR Spectrum , 33 , 885 – 914 . Google Scholar Crossref Search ADS WorldCat Derigs , U. , Pullman , M. & Vogel , U. ( 2013 ) Truck and trailer routing problems, heuristics and computational experience . Comput. Oper. Res. , 40 , 536 – 546 . Google Scholar Crossref Search ADS WorldCat Dror , M. & Trudeau , P. ( 1989 ) Savings by split delivery routing . Transport. Sci. 23 , 141 – 145 . Google Scholar Crossref Search ADS WorldCat Fourer , R. , Gay , D. M. & Kernighan , B. W. ( 2003 ) AMPL. A Modeling Language for Mathematical Programming , 2nd edn. United States of America : Thomson Brooks/Cole . Google Preview WorldCat COPAC Golden , B. L. & Assad , A. A. ( 1986 ) Vehicle routing with time-window constraints: Algorithmic solutions . Amer. J. Math. Management Sci. , 6 , 251 – 260 . WorldCat Ho , S. & Haugland , D. ( 2004 ) A tabu search heuristic for the vehicle routing problem with time windows and split deliveries . Comput. Oper. Res. , 31 , 1947 – 1964 . Google Scholar Crossref Search ADS WorldCat Hoff , A. ( 2006 ) Heuristics for rich vehicle routing problems . PhD Thesis , Molde University College , Norway . Google Preview WorldCat COPAC Hwang , C. L. & Masud , A. S. M. ( 1979 ) Multiple Objective Decision Making, Methods and Applications: a State-of-the-Art Survey . Berlin Heidelberg : Springer . Google Preview WorldCat COPAC Kirkpatrick , S. , Gelatt , C. D. & Vecchi , M. P. ( 1983 ) Optimization by simulated annealing . Science , 220 , 671 – 680 . Google Scholar Crossref Search ADS PubMed WorldCat Kuznietsov , K. A. , Gromov , V. A. & Skorohod , V. A. ( 2016 ) Cluster-based supply chain logistics: a case study of a Ukrainian food distributor . IMA J. Manage. Math ., 28 , 553 – 578 . WorldCat Labbé , M. , Laporte , G. & Mercure , H. ( 1991 ) Capacitated vehicle routing on trees . Oper. Res. , 39 , 616 – 622 . Google Scholar Crossref Search ADS WorldCat Lin , S.-W. , Yu , V. F. & Lu , C.-C. ( 2011 ) A simulated annealing heuristic for the truck and trailer routing problem with time windows . Expert Syst. Appl. , 38 , 15244 – 15252 . Google Scholar Crossref Search ADS WorldCat Mole , R. H. & Jameson , S. R . ( 1976 ) A sequential route-building algorithm employing a generalised savings criterion . Oper. Res. Quart. , 27 , 503 – 511 . Google Scholar Crossref Search ADS WorldCat Mullaseril , P. , Dror , M. & Leung , J. ( 1997 ) Split-delivery routing heuristics in livestock feed distribution . J. Oper. Res. Soc. , 48 , 107 – 116 . Google Scholar Crossref Search ADS WorldCat Nunes De Castro , L. ( 2006 ) Section 3.3: Hill climbing and simulated annealing . Fundamentals of Natural Computing: Basic Concepts, Algorithms, and Applications . Computer Information Science Series . United States of America : Chapman Hall/CRC Press . Osman , I. H. & Laporte , G. ( 1996 ) Metaheuristics: a bibliography . Ann. Oper. Res. , 63 , 511 – 623 . Google Scholar Crossref Search ADS WorldCat Perboli , G. , Gobbato , L. & Maggioni , F. ( 2017 ) A progressive hedging method for the multi-path traveling salesman problem with stochastic travel times . IMA J. Manage. Math. , 28 , 65 – 86 . Google Scholar Crossref Search ADS WorldCat Ruiz , R. , Maroto , C. & Alcaraz , J. ( 2004 ) A decision support system for a real vehicle routing problem . Eur. J. Oper. Res. , 3 , 593 – 606 . Google Scholar Crossref Search ADS WorldCat Tavakkoli-Moghaddam , R. , Safaei , N. , Kah , M. & Rabbani , M. ( 2007 ) A new capacited vehicle routing problem with split service for minimizing fleet cost by simulated annealing . J. Franklin Inst. , 344 , 406 – 425 . Google Scholar Crossref Search ADS WorldCat Yang , S ., Jiang , Y. & Nguyen , T. T . ( 2013 ) Metaheuristics for dynamic combinatorial optimization problems . IMA J. Manage. Math. , 24 , 451 – 480 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Routing problems in agricultural cooperatives: a model for optimization of transport vehicle logistics JF - IMA Journal of Management Mathematics DO - 10.1093/imaman/dpy010 DA - 2019-09-02 UR - https://www.deepdyve.com/lp/oxford-university-press/routing-problems-in-agricultural-cooperatives-a-model-for-optimization-N0RxAzEwNN SP - 387 VL - 30 IS - 4 DP - DeepDyve ER -