#### Abstract

This paper addresses a biobjective production-distribution planning problem. The problem is formulated as a mixed integer programming problem with two objectives. The objectives are to minimize the total costs and to balance the total workload of the supply chain, which consist of plants and depots, considering that it represents a company vertically integrated. In order to solve the model, we propose an adapted biobjective GRASP to obtain an approximation of the Pareto front. To evaluate the performance of the proposed algorithm, numerical experimentations are conducted over a set of instances used for similar problems. Results indicate that the proposed GRASP obtains a relatively small number of nondominated solutions for each tested instance in very short computational time. The approximated Pareto fronts are discontinuous and nonconvex. Moreover, the solutions clearly show the compromise between both objective functions.

#### 1. Introduction

Nowadays, supply chains and manufacturing systems must be responsive, agile, and flexible to adapt to changes and requirements of the markets. Globalization and the fierce competition in the markets impose more challenges to such systems, increasing their complexity. The emergence of disruptive technologies in which data can be captured by sensors in real-time represents an opportunity to improve decision-making, but it is necessary to integrate it with proper analytical models to support planning and scheduling decisions along the supply chain and manufacturing systems.

One of the main goals of Supply Chain Management is the coordination of the different echelons and the related planning decisions such as supply, manufacturing, inventory management, and distribution [1]. A typical approach is to address previous planning decisions as independent problems and obtain an overall solution hierarchizing the order of the decisions and solving the problems sequentially. This is motivated by the complexity of modeling them into a single integrated problem, but it has some drawbacks as it does not consider the interrelations among the decisions. Hence, the use of integrated problems provides a more effective and efficient way to obtain solutions for the underlying real problem and reduce total landed costs.

There is a very close linkage between production and distribution planning decisions [2], reason for which these problems have been addressed in the literature as a joint production-distribution planning problem. Several contributions can be found in the literature addressing coordination in the supply chain in which production and distribution problems are integrated. In this paper, we address a production-distribution planning problem for a single product faced by a manufacturing company. According to the specific characteristics of the company, the problem is formulated as a biobjective mixed integer programming model. The objectives consider the minimization of the total costs incurred by the company and the balancing of the total workload assigned to the plants and depots of the company. This is justified as we are considering the case of a vertically integrated supply chain. Hence, the company is interested not only on the minimization of total costs, but also on a better resource utilization. Balancing the workload content in the supply chain contributes to achieve a better utilization of the resources already installed in the supply chain network, in addition to avoid congestion at the plants and depots. To solve the model, we propose an adapted biobjective GRASP to obtain an approximation of the Pareto front.

##### 1.1. Literature Review

One of the earliest contributions on production-distribution planning in the literature is provided by [3]. They address a multicommodity, single-period production-distribution problem and propose a solution methodology based on Benders decomposition. Another early contribution is provided by [4]. They propose a methodology that considers the relationship between a multistage production and distribution system, capturing the stochastic and dynamic interactions. The proposed model is composed by a material control submodel, a production submodel, a stockpile inventory submodel, and a distribution submodel.

There exist several survey articles in the literature that provide a framework or taxonomy of related production-distribution planning problems. We briefly review the existing surveys and cite some of the contributions found in the literature, with special emphasis on the most recent ones found in the literature. An overview of the literature related to the coordination of production planning among multiple plants assuming vertically integrated firms is provided by [5]. A review of the literature related to coordinated planning between two or more decisions/echelons of the supply chain is presented by [6]. One of the categories analyzed correspond to production-distribution planning problems. They propose a taxonomy of the contributions as buyer-vendor coordination, production-distribution coordination, and inventory-distribution coordination. A taxonomy of the major decisions for integrated production-distribution planning is provided by [7]. They emphasize the complexity of the operations and make a reference to the emergence of philosophies as a response to competitive pressures of markets, such as Total Quality Management and Flexible Manufacturing Systems and Supply Chain Management. In this case, the aim is to achieve flawless quality of products and quick response and agile manufacturing at low cost and an efficient distribution of products to the consumer points with low stocks.

In addition, a review of the strategic production-distribution models is presented by [8], emphasizing global logistics systems and mixed integer programming models. They classify the contributions into four categories: previous reviews, optimization models, additional issues for modeling, and case studies. Another review of integrated analysis of production-distribution systems is presented by [9]. The aim of their review is to determine how the different logistics aspects can be included in an integrated analysis and what competitive advantages can be obtained when integrating production and distribution decisions. They include only those contributions that explicitly address transportation systems and provide a classification of distribution-inventory models pointing out the contributions found in the literature [10–15].

Another related survey is presented by [16]. They categorize the contributions and provide a framework for the existing literature and define a research agenda. A classification of production-distribution scheduling is proposed in [17]. They point out that given that production and distribution stages in a supply chain are often linked by an inventory stage, contributions in the literature often involve inventory decisions and noticed that the literature that consider scheduling decisions has recently attracted the attention of researchers in the literature [18–21]. They divide existing models into five classes: (i) models considering individual and immediate delivery, (ii) models with batch delivery to a single customer by direct shipping method, (iii) models with batch delivery to multiple customers by direct shipping method, (iv) models with batch delivery to multiple customers by routing method, and (v) models with fixed delivery departure dates [22–26].

On the other hand, a review the existent approaches based on the modeling approach and the solution techniques employed is presented by [27]. They provide a taxonomy of the contributions in the field. The most recent review is provided by [28] that considers the evolution of production-distribution planning problems towards the collaborative production-distribution planning in Supply Chain Management. They analyze how the integrated approaches have an impact on the industry and their real requirements.

In terms of the solution techniques, it is possible to find mathematical programming-based approaches or solved directly by a commercial solver (e.g., [3–5, 10–16, 21, 23, 29–37]); a vast number of contributions that propose heuristic techniques, motivated by the complexity of the resulting problem [18, 19, 22, 26, 32, 38–47]; and simulation-based approaches to tackle the problem [48–51].

Another important feature is related to the objectives that are optimized. Several contributions consider a single-objective optimization approach [2, 22, 29–33, 39, 41–46, 49, 52, 53]. However, more recent contributions are capturing more realistic features of planning decisions by the incorporation of multiple objectives to be simultaneously optimized [40, 51, 54]. The typical objectives considered are profit maximization, costs minimization, and service level offered to customers. For instance, in [54], the minimization of total costs and the minimization of total tardiness of products transported to the distribution centers are considered as objectives. In [55] the focus of the proposed model is on customer’s service level, measured as the mean delivery time and total transportation cost. On the other hand, in [56], minimizing total operation cost of the vehicles to deliver the products and the waiting time of the customers is considered as an objective; hence they are more focused on the distribution problem, although they also consider production scheduling decisions.

None of previous contributions have incorporated workload balance as an objective to optimize. Workload balancing has been considered in decision planning problems related to vehicle routing problems [57–60] and machine scheduling problems [61]. Balancing the workload content of the supply chain is an important element when considering a vertical integrated company. This is meaningful as the company aims to have a good utilization of its assets. Also, having a plant with overutilization generates congestion while others may be underutilized. To the best of the knowledge of the authors, this is the first paper in the literature that incorporates workload balance as an objective of the integrated production-distribution planning problem.

Another contribution of this work is related to the specific features of the solution approach. Only few papers have adapted GRASP to solve biobjective problems. Furthermore, we propose a clear integration of both objective functions when constructing solutions. Moreover, in the construction method, we hybridize mathematical programming techniques with heuristic ones. The local search method explores typical neighborhoods for routing problems but also mathematical programming techniques are used to complete a movement.

This article is organized as follows: Section 2 presents the proposed mathematical model. Section 3 describes the heuristic algorithm and Section 4 presents the numerical results showing the suitability of the proposed algorithm to approximate the Pareto front. Finally, Section 5 presents the final remarks and recommendations for further research areas.

#### 2. Biobjective Optimization Model

The problem considered in this manuscript is related to production and distribution planning decisions of a single product faced by a manufacturing company. For the production part, the company is responsible for manufacturing products in plants, respecting the capacity of the plants. Also, the company must satisfy the demand of all customers that are served by a depot in which products are consolidated, similar as in a standard transportation problem. On the other hand, for the distribution part, the company has to determine the routes to serve the demand of customers that are assigned to each depot. For this, the classical assumptions considered in vehicle routing problems are considered. For instance, the routes should begin and finish in the depot; there is a maximum duration time established by each route, a customer is visited once by a single vehicle and a single vehicle is associated with a route, and the demand of all customers is satisfied. Furthermore, a homogeneous fleet of vehicles is considered in terms of capacity. In terms of network design, it has been established that depots are not connected to other depots and that a depot can have multiple vehicles assigned, each serving a single route. It is also assumed that depots are uncapacitated, while plants have capacity constraints.

The planning decision of the company aims to optimize simultaneously two objectives while respecting previously described constraints. The first objective is to minimize total costs (production and distribution), and the second objective is to balance total workload of plants and depots. Balancing the workload content at each facility contributes to achieve better resource utilization over the network; and, accordingly, it is expected to achieve less congestion at the facilities such as plants and depots. Under this context, a straightforward solution would be to equally split the production among plants and customer’s demand among depots. However, this solution would not necessarily minimize production and distribution costs. Hence, the company is facing a production-distribution problem that considers typical constraints of a transportation problem merged with a multidepot vehicle routing problem (MDVRP).

Prior to describing the mathematical model of the proposed problem, the sets, parameters, and decision variables are defined. Here consider as the set of plants, the set of depots, the set of customers, and the set of vehicles. Let be the set of the arcs connecting depots and customers. The customers served by depot are included in and customers included in vehicle are included in . Also, let be the set of vehicles used in depot .

Note that there are distinct costs involved in the model. The operating cost associated with manufacturing a product in plant and ship it to depot is denoted by . The distribution cost associated with the shipment of a product from depot or customer to depot or customer , where is denoted by . Each customer has a demand , and the manufacturing capacity at each plant is . Furthermore, the time required for travelling through an arc is denoted by , the required time to unload customers demand , from any vehicle is given by , the capacity of vehicle is , and the maximum duration time of vehicle in a route is .

In this problem, are binary decision variables and are continuous ones. The fact that implies that vehicle uses arc. Also, represents the amount of products manufactured in plant that are shipped to depot .

To balance workload in plants and depots, the deviations with respect to the ideal case are computed. The latter refers to the case in which the workload is equally distributed in plants and depots, respectively. For the plants, the ideal case is when all of them are at the same capacity level; that is, the production at each plant is proportional to their capacity such that total production satisfies the total demand of customers. For instance, plants with greater capacity may manufacture more products than plants with lower capacity. For the depots, the ideal case is when the customers demand is equally distributed among all the depots. Notice that, in this case, the balancing of the workload at the plants and depots is modeled based on a deviation value of the products manufactured at each plant and consolidated at each depot, respectively. This resembles the objective of the company to have a target value of workload at each facility.

Based on these ideas, we compute the maximum deviation for the plants as follows:The maximum deviation associated with a depot is computed as follows:Note that the quantity of products manufactured in plants is greater than the amount of products consolidated in depots. Thus, the value given in (1) may be greater than the value given in (2). In this case, we are balancing the workload content only at the plants nor in the depots. To overcome this issue, an value is included to standardize the units in the sum.

The resulting mixed integer biobjective programming problem is defined as follows: Equation (3) represents the first objective function, which aims to minimize the sum of the operating cost (manufacturing and shipping) and the cost from shipping products from depots to customers. The second objective function is expressed by (4), which aims to balance the workload in depots and plants. It expresses the corrected sum of the maximum deviation from the ideal case split. The second objective function aims to find a solution with small deviations, that is, with a balanced workload content in plants and depots. Production constraints are expressed in (5) and (6). Constraints (5) ensure that the amount produced in each plant does not exceed its maximum manufacturing capacity. Constraints (6) guarantee that the demand of customers assigned to each depot is satisfied. Distribution constraints are expressed in (7)–(15). Visiting customers once by a single vehicle is forced in constraints (7) and (8). Constraints (9) and (10) determine the assignment of customers to depots. Despite the fact that, in order to balance workload all depots must be used, to minimize costs it could be desired to use only a subset of depots. Constraints (11) and (12) prevent the interdepot flow. This is included as shipments between depots are not allowed. Subtours are prohibited in constraints (13); constraints (14) determine the vehicle’s capacity constraints, and the maximum duration of a route performed by a vehicle is defined by constraints (15). Equation (16) expresses the binary nature of the routing variables, and (17) establishes the nonnegativity of the manufacturing variables.

#### 3. Description of the Proposed Biobjective GRASP

GRASP is a metaheuristic that has been proven to be efficient when solving complex combinatorial optimization problems [62]. It consists of two main methods. The first method is devoted to construct an initial solution and the second method concerns to the local search, which aims to improve the current solution. As GRASP is a multistart metaheuristic due to the randomized procedure that is implemented in the construction method, it is possible to construct multiple solutions that are further improved during the local search method. Hence, the two methods are iteratively repeated until a stopping criterion is reached. The best solution in terms of objective function value is returned as the output of the algorithm.

In addition to the successful use of GRASP to solve single-objective problems, it has also been successfully applied to approximate the efficient Pareto front of different multiobjective optimization problems [63–68]. The previous cited references adapted specific components of GRASP to the corresponding problems under study. However, a taxonomy to standardize the adaptation of GRASP to multiobjective problems is provided by [69]. Thus, we follow previously referred taxonomy during the description of the proposed algorithm. The main elements of GRASP for biobjective problems are the following: a greedy function is considered when constructing a solution, but, in case of having multiple objectives, a criterion must be selected to identify the guiding function, the local search method considers the multiple objectives to evaluate if the current solution has been improved or not, and a set of efficient solutions is provided when the algorithm reaches some stopping criterion instead of returning the incumbent solution.

In the proposed biobjective GRASP a pure-random construction is done; that is, one of the two objective functions is randomly selected and all the solution’s construction is guided by the greedy function associated with that objective function. Furthermore, a pure strategy is considered in the local search with respect to one objective function randomly selected (note that it may be different to the objective function considered in the construction). In this case, an improved solution may have better objective function value for the considered objective function but it may be worsened in the other objective function value. One important feature of our proposed algorithm is that the objective function considered in the local search may be different than the objective function used when constructing the solution. Usually, in pure strategies, the objective function considered during the construction and local search is the same.

##### 3.1. Construction Method

The first step of the construction method is to randomly select one objective function. Here, both objective functions are equally likely to be selected; that is, if , then the objective function that concerns to total costs is selected (see (3)), and in case when , the objective function related to the workload balance is considered (see (4)), where is a random number between 0 and 1. For any of the two objective functions, the construction method can be divided to three phase: clustering, routing, and manufacturing.

###### 3.1.1. Clustering Phase

In the clustering phase, the allocation of customers to depots is made by following random-greedy strategies, in which an allocation problem is solved depending on the objective function under consideration. If (3) is being considered, then the auxiliary problem defined by (18)–(20) is optimally solved by CPLEX, where is a random number.Since, we are concerned in constructing a good solution in terms of the total cost, the minimization of the distribution costs from customers to depots is pursued in (18), considering that customers must be allocated only to a single depot. A perturbation in the distribution costs is considered to allow certain degree of randomness in the clustering phase.

In the other case, when (4) has been selected, the allocation of customers to depots such as a balanced workload is obtained. In this case, the problem defined by (21)–(25) is optimally solved by CPLEX, where is a predefined constant value and is a random number.Bearing in mind that we are seeking to construct solutions with balanced workload, constraints (22) and (23) will act as balancing constraints but considering a tolerance value . A small value of implies that more balanced clusters are given. Analogously, greater values of allow the depots to have less consolidated customers demand associated with them. Note that if , the ideal case occurs (i.e., when all depots have the same workload), but this case may be infeasible due to the unique allocation of customers to depots imposed by (24).

###### 3.1.2. Routing Phase

Once customers have been allocated to the depots, it is necessary to design the routes to distribute the products to the customers. For this, greedy-randomized adaptive strategies are considered. To create the routes, we consider that a single vehicle departs from a depot and the route is created by incorporating the customers to be visited, assigning one customer to each depot at a time, that is, in a parallel manner. Regarding the available number of vehicles in each depot and their capacities, the same assumptions of [47] are considered.

Note that the routing is made for customers allocated in the same depot. Hence, it is evident that the workload in each depot holds. Based on that, we create routes with low distribution cost as it is not necessary to take into account the other objective function. Thus, for any of both objective functions, the routes are created in the same manner.

To add a customer in a route, the associated cost in the guiding function that summarizes its contribution is computed. As guiding function, we consider the cost of adding customer after customer in the current route, in which the depot is denoted by ; that is, . The latter is computed for all customers in a particular depot . Then, considering the same depot, a Restricted Candidate List () is created. One criterion considered to add a customer into the RCL is when , in which and denote the minimum and maximum cost of adding a customer. Also, is a parameter that regulates the degree of greediness considered in . Note that when , customers are selected in a greedy manner, and if customers are randomly selected. Other two important criteria that customers must meet to be included in are that vehicle’s capacity and maximum duration of the route may not be exceeded. When the are formed (for all ), a customer is randomly selected to be added in the current route of depot . As it is mentioned above, one customer per depot is added at a time, that is, a parallel routing. This strategy continues until all customers are in a route. It is worth clarifying that if there are no customers that satisfy all the criteria mentioned above, then . In this case, a new vehicle is used to begin a route.

###### 3.1.3. Manufacturing Phase

The last phase of the construction method corresponds to the manufacturing decisions in plants. The decision taken in this phase is the amount of products manufactured in each plant. This resulting problem may be seen as a classical transportation problem, in which different structures are considered depending on the current target objective function. For example, if the first objective function is being considered, it is aimed to manufacture the products in plants with lower operating costs. On the other hand, if the second objective function is under consideration, it is aimed to balance workload of the plants according to their production capacity without regarding the operating cost.

Therefore, if (3) is considered, the following auxiliary transportation problem is optimally solved by CPLEX.The problem defined by (26)–(29) attempts to minimize operation costs in the plants, while respecting the production capacity in the plants (see (27)) and satisfying the consolidated demand in each depot (see (28)).

If (4) is selected as the objective function of interest, the auxiliary problem must be adapted in such a way that the workload is balanced. The resulting problem is presented next, where is a predefined constant value. The problem defined by (30)–(33) is optimally solved by CPLEX. Note that (31) equally distributes the production among the plants.

After one solution has been constructed, the corresponding values for both objective functions are computed. That is, if the solution was constructed based on the minimization of costs, the corresponding value associated with the workload balance is evaluated, and vice versa.

##### 3.2. Local Search Method

The implemented local search method attempts to improve the constructed solution in terms of one of both objective functions; that is, one objective function is randomly selected. The local search is mainly performed in the routing phase with neighborhoods named remove/insert and interchange, and , respectively. Both neighborhoods can consider intradepot and interdepots movements. During the exploration, the best improvement criterion is assumed. The considered movements are adapted from the ones described in [47], but they are further explained for each objective function.

If the first objective function is being used as a guiding function (see (3)), then and are sequentially explored. First, is considered, in which a predefined small number () of customers with higher costs are removed from each route. Then, one of those customers is randomly selected to be inserted in the best position of any route that do not violate any constraint (capacity of the vehicle or maximum duration of the route). A particular customer can be inserted into a route of the same depot or into a route of a different depot. In the case that the insertion procedure is performed in a route belonging to a different depot, the consolidated demand of the involved depots changes. Thus, the manufacturing decision is updated by solving the transportation problem that minimizes (26) subject to constraints (27)–(29). That problem aims to minimize the operating costs associated with the manufacturing process without being aware of the workload balance in the plants. On the other case, if the consolidated demand in the depots remains the same, there is no need to updating the manufacturing decisions. The exploration of terminates when all the customers are inserted into a route.

After that, is explored. Here, two customers from different routes are interchanged with each other. An interchanged customer occupies exactly the same position of the other interchanged customer, but in the corresponding route. As before, an interchange may or may not modify the consolidated demand in the involved depots. If the consolidated demand is modified, then the problem defined by (26)–(29) is optimally solved. In the other case the manufacturing phase does not need to be updated.

Now, the case when the second objective function (see (4)) is considered is described. The guiding function to improve the current solution is (2). For , the customers removed from each route are randomly selected. Then, those customers will be inserted into a route (possibly from a different depot) that improves the workload balance of customers demand in the depots. Note that there is no significant difference regarding the route of the same depot in which a customer is inserted. Thus, routes associated with the same depot are randomly selected. In this neighborhood, the manufacturing decision must be updated after the best improvement has been conducted. The problem defined by (30)–(33) is optimally solved to achieve the latter. The exploration of is adapted in an analogous manner compared to the one described for the other objective function. That is, customers are interchanged between routes but focusing in the workload balance objective. By doing this, we can evaluate the convenience of the movement performed. After the best interchange has been conducted, the problem defined by (30)–(33) associated with the manufacturing phase is solved.

Without regarding the objective function considered to conduct the exploration of and , we measure the improvement in the solution by checking if it is included in a nondominated set (). An improved solution is included in , if there is no solution in that simultaneously outperforms both objective functions of the current solution. The stopping criterion is a predefined number of iterations without updating the set.

A diagram flow describing the steps of the algorithm is shown in Figure 1.

#### 4. Computational Experimentation and Results

Computational experimentation is performed to measure the performance of the proposed biobjective GRASP algorithm. A set of ten instances previously used in the literature to solve a similar production-distribution planning problem is considered. The instances used in [70] were adapted to create instances for the problem herein studied. To adapt the instances, acquisition costs are omitted and the value that standardize the magnitude of workload in depots and plants is included. As indicated in Section 3.1.1, the values for and were randomly defined in a range of and . The information regarding the number of customers, depots, and plants is shown in Table 1.

The proposed algorithm was implemented in Microsoft Visual Studio 2012 coded in C++. The version of CPLEX used to solve the auxiliary problems is 12.6.1. All the computational experimentation was performed in a personal computer Intel (R) Core processor with a speed of 2.60 GHz and 4 GB of RAM. Since the biobjective GRASP has stochasticity in some parts, each instance was solved five times. Table 2 summarizes the results obtained from all the runs. The first column indicates the instance and the second column shows the average of the number of nondominated solutions in the five runs. The third and fourth columns display the minimum and maximum number of nondominated solutions in from all the runs. The last column shows the algorithm’s average required time (in seconds) for solving each instance.

From the results shown in Table 2, it can be seen that a relatively small number of nondominated solutions are found by the proposed algorithm (from (9) to (28)). Particularly, for Instance , a small number of nondominated solutions was found in each run. On the contrary, the largest number of nondominated solutions was obtained, for instance, . Also, in five out of the ten instances, the minimum and maximum number of nondominated solutions are very similar, which indicates that the algorithm performs steady (see instances , , , , and ). Despite the relative small number of nondominated solutions obtained, the latter fact validates that the proposed algorithm has a good performance. In addition, note that this combinatorial problem tends to have a relative small number of solutions in the Pareto front due to the structure of the problem.

It is convenient to highlight that the computational times are very short for a complex operational problem of this nature. Hence, the stopping criterion may be changed to explore more solutions with the aim of increasing the number of nondominated solutions.

To illustrate the approximated Pareto front of the tested instances, the nondominated solutions are plotted in Figures 2–11. The run with more solutions in the nondominated set is plotted for each instance. In these figures, the horizontal axis represents the objective function associated with the minimization of total costs (see (3)) and the vertical axis corresponds to the workload balance (see (4)). Note that a solution with low cost tends to be less balanced (has higher deviation in the workload); on the contrary, a solution highly balanced is associated with higher costs.

It can be seen from Figures 2–11 that the obtained approximations of the Pareto front are not covering all the space as expected. Some regions without nondominated solutions are observed. As reported by [71], depending on the modeling approach for balancing workload, it is possible to find less number of solutions for the approximated Pareto front. They tested four different schemes to model the workload balance objective function for a capacitated vehicle routing problem and results show that the number of Pareto optimal solutions varies substantially. A less number of Pareto solutions is found when the balancing constraint is modeled as the total difference of travel distances of a route with respect to a desired target or average value. Note that this is the approach that we are employing to model the deviations of the workload assigned to the plants and depots.

For example, for Instances , , , and , the nondominated solutions are very polarized. In other words, there are many nondominated solutions regarding the first objective and less solutions that favor the workload balance. But the tradeoff between both objectives triggers a discontinuity in the Pareto front. That is, for a small increment in the value of the first objective function (axis ), the value associated with the second objective function (axis ) significantly decreases close to zero. Thus, the Pareto front is discontinuous, which is common in multiobjective combinatorial problems [72]. Also, for almost all the instances the approximation of the Pareto front is nonconvex, which emphasizes the difficulty of the problem herein studied.

#### 5. Conclusions

In this paper, we introduce a biobjective production-distribution planning problem for a two echelons vertically integrated supply chain. The model considers both the minimization of total costs and the workload balance among plants and depots involved in the supply chain. The latter objective has been commonly used in vehicle routing problems and manufacturing systems. Due to the vertical integrated nature of the company, balancing the workload among the facilities is an important performance indicator to achieve. Thus, the introduced model contributes to the literature by simultaneously considering these two objective functions in a production-distribution planning problem.

Because of the complexity of the problem herein studied, we propose a heuristic approach based on the well-known metaheuristic GRASP to obtain an approximation of the Pareto front. GRASP was adapted to consider both objective functions during the construction and local search methods. To evaluate the performance of the proposed algorithm, we employ a set of instances from the literature used for a similar problem. Nevertheless, the results reported in the literature that correspond to the used instances cannot be compared with our obtained results due to the structural differences among the respective considered models. Results indicate that the Pareto front is discontinuous for almost all the tested instances. According to [72], this fact is common in combinatorial biobjective problems with constraints that cause discontinuous solution space. In the cases when the approximation of the Pareto front are discontinuous, the nondominated solutions are polarized in the extreme of the regions; that is, the solutions are very good for one objective function but very bad for the other objective function (this can be seen from Figures 2–11). In addition, the obtained approximation of the Pareto front is nonconvex for almost all the tested instances and, as remarked by [71], the number of nondominated solutions obtained depends on the manner that the workload balance objective is modeled. In this paper, we modeled it aiming to have a workload with less deviation from a target value (the one associated with the equally balanced case). Hence, the Pareto front does not cover all the solutions space and is discontinuous and nonconvex due to the structure of the considered model.

As future research, we propose to apply path relinking to the obtained nondominated solutions to validate the discontinuity of the approximated Pareto fronts. If no significant increase in the number of nondominated solutions is obtained, we may experimentally conclude the discontinuity. Moreover, different modeling approaches to balance the workload can be used (as in [71]). By doing the latter, we may theoretically evaluate their appropriateness and applicability to the production-distribution problem considered in this paper.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The research of the first two authors has been partially supported by the Mexican National Council for Science and Technology (CONACYT) through Grant SEP-CONACYTCB-2014-01-240814. Also, the second author acknowledges the program of Professional Development of Professors with Grant PRODEP/511-6/17/7425 for research stays during his sabbatical year.