Two Efficient Meta-Heuristic Algorithms for the Robust Inventory Routing Problem with Backhaul

: The inventory routing problem (IRP) involves the integration and coordination of two components of the logistics value chain: inventory management and vehicle routing. Therefore, consideration of this issue can be effective in decision making of the organi zation and will lead to lower costs or other goals. Our objective in this article is to examine a new inventory-routing model and solve it with meta-heuristic methods. For more flexibility of the model, and approaching the real world, the model of this article is considered multi-period and multi- product. Also, two objective functions, including minimizing system costs and transportation risk, are included in this model . Given that the main parameter of the model, that is, demand, is uncertain, we have used a robust optimization approach to solve it, and since this model is in the classification of NP-Hard problems, we have used two meta-heuristic algorithms consisting of non-dominated sorting genetic algorithm (NSGA-II) and a multi- objective imperialist com petitive algorithm (MOICA) . By examining the model in two deterministic and robust conditions, according to two criteria, the mean values of the object ive function and its standard deviation, it has been determined that in almost all cases, the robust optimization model produces better solutions. Also, between the two meta-heuristics method, the NSGA-II algorithm has shown better quality according to the mentioned criteria. Obviously, taking into account the different features of a model increases its efficiency. But this, obviously, makes the model even more complex. However, this complexity of models can work like a real system. Our attention in this article has been to this subject. To analyze such models, exact methods do not have the required performance and paying attention to heuristic and meta-heuristic methods is very effective. In this paper, a robust optimization and meta-heurictic approaches focus on these goals.


INTRODUCTION
In recent years, the increasing development of the chains and the created competition among them, on the one hand, and the development of the information management and the greater awareness of the companies of their chain performance, on the other hand, have led to a great consideration for the coordination, co-operation, and integration of the various elements of the supply chain in order to achieve the competitive advantage. The literature, in this area, shows that coordination in the production, inventory, and transportation management will reduce the costs and increase the total service level of the chain. In the supply chain, the different products are often produced and consumed in different locations. Hence, transportation is one of the key elements in the supply chain and it may generate one of the major parts of the costs in the supply chain. The design of the transportation network affects the performance of the supply chain. The proper design of the transportation network in a supply chain can result in spending a low cost to achieve a satisfactory level of accountability [1].
In a supply chain, different goals are considered. One of the most important of these goals is the cost debate. In connection with this, transportation plays an important role. That is why designing a proper transportation network is very important in achieving the goals of the chain. According to this discussion, vehicle routing problem (VRP) was formed as an optimization problem. In VRP, we are looking for ways to transfer goods so that they have the lowest costs at any given period. Over time, the companies realized that merely considering a routing segment would not reduce costs in the long-term planning, and the attention to the logistics segment and inventories could also add to the optimality of the problem [2].
The vender-managed inventory (VMI) system is one of the most recent examples in creating added value through logistics services. This policy is often described as a win-win situation for both, the supplier and the customer, leaving customers free from high investment on the inventory and the complexity of the inventory control. Therefore, the VRP gains more significance when it also includes the simultaneous management of the inventory. In this case, the raised issue would be the inventory routing problem (IRP). In the IRP issues, the planning horizon is usually over a period, and the goals of the chain are examined throughout the horizon. That is, we are looking to reduce costs (or other objectives) at the end of the horizon [3].
In the case of VRP, customers send their order to the seller, and the vendor assigns these orders to the vehicleroute, so that the total distance traveled will be minimized while satisfying the orders. But in the IRP issue, the vendor decides on the amount of product that should be sent to each customer daily, not the customer himself. In other words, every day from the planning horizon, the vendor decides which customers to receive and how much they receive.
Therefore, the inventory-routing problem generally involves three choices: • When will each customer be served? • How much product will be sent to customers who are served to them? • What routes should be considered for sending?
Numerous studies and analyses have been carried out on IRP previously; for example comprehensive review has been presented in [4].
IRP issues have different varieties in which the inventory-routing problem with backhauls (IRPB) has been considered. At IRPB, the customers are divided into two groups of linehauls (delivery of the goods) and backhauls (receiving the goods). The linehauls customer needs a number of goods to be delivered, and the backhaul customer also needs a number of goods to be received. This issue has not yet been addressed, in the literature, of the IRP.
In further research, we study the literature of review. Attracted articles include articles that have been presented in valid journals and conferences in the field of IRP between 2000 and 2019. In Tab. 1, a summary of the review articles of this research is presented according to our classification of these articles on the number of objective functions (single-objective or multi-objective), the nature of model data (certain or uncertain), the number of products (single or multi), the planning horizon (single or multiple) Periodic) and a brief explanation of the method of solving the model is given. Considering the parameters of each model in real terms, it can have a lot of effects on its performance. One of those parameters in IRP issues is demand that we assumed to be uncertain. To solve such uncertain models, there are various methods, such as stochastic programming, dynamic programming, fuzzy methods, and robust optimization. But given the fact that the accurate distribution of the parameter of demand is not clear, and if we can guess the range of the change of this parameter, the best approach is robust optimization, which tries to preserve as much as possible the optimality under uncertain parameters.
Robust models can be classified based on conditions and attributes of models. In Tab. 2, this topic is mentioned. It should be noted that scenario-based models are appropriate when uncertainty occurs in a discrete or scenario form. There is the possibility of combining robust optimization models with other approaches to dealing with uncertainty, such as the fuzzy approach. Table 2 Classification based on conditions of models Approaches Some related research Scenario-based stochastic robust programming [19] Robust programming based on closed convex uncertainty sets [20] Fuzzy robust programming [21] This paper examines a multi-product and multi-period IRPB with the aim of minimizing holding and transportation costs with risk. The IRP is studied in a supply chain consisting of two parts -a distributer, and a set of retailers with direct delivery strategy. The demand for this model is non-deterministic, and robust optimization approach has been developed. Also, the metaheuristic methods, namely NSGA-II and MOICA, have been used to solve it. Also, considering the multi-objective IRPB model in this paper, it has not yet been seen in the IRP literature and therefore the research innovation is summarized in this case.
In this paper, the following sections are presented: model description in section 2, robust optimization in section 3, the solution approaches in section 4, and discussion in the final section.

MODEL DESCRIPTION
In this section, the mathematical model of the research paper is described. The structure of this model is a mix integer programming (MIP). It should be noted that Arab et al. [22] presented the basic model of this research in year 2020. The following assumptions are also considered: • Model is single-depot, which has to supply all demands of all customers. • Distribution fleet is not homogenous. • Distances between customers are specified. • Inventory shortage is not allowed.
In the following, before describing the mathematical model, the parameters and variables of the model are described. . .

Symbols
, 1; , , , The first objective function covers system costs, fixed routing costs included, transportation and maintenance cost. Also, the second objective function minimizes the risk between the nodes. Eq. (3) and Eq. (4) express the inventory balance for customers on a round trip with respect to their demand, respectively. Eq. (5) and Eq. (6) represent the difference between the input and output of each node for customers on a round trip. Eq. (7) states that each customer is visited by a vehicle maximum once during each period. Eq. (8) and Eq. (9) indicate that each vehicle starts at the central depot and returns to that after each trip. Eq. (10) shows the continuity of the travel path. Eq. (11) indicates that customers on the linehaul trip are to be visited and provided with service before customers on the backhaul trip. Eq. (12) indicates the customer's inventory capacity. Eq. (13) is used to sub-tour elimination from vehicle routing problems, and Eq. (14) and Eq. (15) show the maximum and minimum load of variable products for each vehicle during a linehaul or backhaul trip. Also other constraints show type of variables other constraints, present assumptions and values of variables.

ROBUST OPTIMIZATION
For the first time, Bertsimas & Sim [20] offered robust optimization for discrete problems. This method is applicable to the linear optimization problem; in the case of uncertainty coefficients exist both in the objective function and in the constraints. This approach also applies to linear continuous optimization problems. We consider the following optimization problem in general:

Model of Data Uncertainty U
Uncertainty parameters can include the coefficients of the objective functions or the coefficients of the constraints. They are referred to below: i. Uncertainty for matrix A: Each of the constraint coefficients a ij , j ∈ N = {1, 2, …, n} is modeled as a non-dependent random variable, with uniform

Robust MIP formulation
Consider the i-th constraint problem as is defined as a set of uncertain coefficients in row i. For each row we define the parameter Γ i that Γ i ∈ [0, |J i |], (means |J i | is the total number of parameters that have changed in the i-th row). In fact, Γ i has the role of adjustment of stability model. Authors showed that, with little probability, all coefficients could be uncertain. We therefore assume that the maximum of [Γ i ] is the number of these coefficients allowed to change, and a coefficient a ij can also change to  Also, due to the symmetric distribution of variables, even if some coefficients that change are also increased from [Γ i ], the optimal solution will be justified with a high probability. Therefore, we call Γ i the level of protection for the i-th constraint.
Also the parameter Γ 0 adjusts the level of robustness in the objective function. Let J 0 = {j|d j > 0}. If Γ 0 = 0, we ignore the impact of the cost deviations, while if Γ 0 = |J 0 |, we are examining all possible deviations, which is actually the most conservative. Therefore, the proposed robust counterpart of Eq. (17) is as follows: If we are interested to make the above model into the linear mathematical model, then the theorem (1) is required. Theorem 1. Given the vector given x * , the protection function of i-th constraint is obtained from Eq. (19).
Which is equivalent to: The proof of the aforementioned theorem is presented in the paper by [20]. By placing the dual of the Eq. (20) in the initial robust model, it can be rewritten as follows: It should be noted that the variables added in the dual robust model (z i , y j , r ij , z 0 , r 0j ) are used to adjust the robustness of solution and apply the protection levels in the model in which r and z are vectors of dual variables in the constraints and objective function that are used for linearization of nonlinear relations.

Formulation of our Robust Model
In our model, we assume that demand is not deterministic. This topic is shown in Eq. (3) and Eq. (4). To counteract this uncertainty, we use Bertsimas robust optimization approach. Therefore p it d and ' p it d appear in a non-deterministic form in the model. Therefore, in Eq. (3) we define the coefficient , and is a number between zero and the number of non-deterministic parameters associated with this row. Or, in other words, the number of changes in demand for product p from the linehaul customers during the planning period. Also in Eq. (4), for backhaul customers, we define such relationships as above.
Since there are no parameters d and d' in the objective function that has uncertainty, we only formulate demand constraints based on the Bertsimas approach. We have a general dual form of it: that: Also: As a result, the final linear robust form will be as follows (Eq. (1') to Eq. (19')):

SOLUTION APPROACHES
Given that the model of this research is bi-objective and the objectives are inconsistent, multi-objective approaches (MOP) may be the best approach to solve it. On the other hand, because the research model is NP-hard, exact solution methods cannot be efficient in a reasonable time. In this case, meta-heuristics methods are very effective. In this paper we use two methods, namely Nondominated Sorting Genetic Algorithm, version2 (NSGA-II) and Multi Objective Imperialist Competitive Algorithm (MOICA). These methods are described below.

NSGA-II
Srinivas & Deb [23] proposed a non-dominated sorting genetic algorithm (NSGA). This algorithm has two main operators: • fast non-dominated sorting • crowding distance. With the help of these operators, this algorithm can be implemented for multi-objective problem, and by creating the Pareto Front, it solves the model. Therefore, the second version of the NSGA is very useful in solving multi objective problems.
The efficiency of the NSGA-II algorithm has been shown in many articles. Therefore, we also use this algorithm and, of course, MOICA algorithm to solve the IRP model in this paper. Also, the pseudo code of the algorithm used in this paper is shown in Fig. 1.
The main points are considered in relation to this algorithm: • For solution representation, we use the matrix with one row and U + W + V − 1 columns, where U and W represent the number of customers and V represents the number of vehicles. • The crossover operator uses the clever single point method and the mutation operator uses the swap method. • Stopping criterion for the proposed algorithm is equal to the maximum number of repetitions of it.
• The selection process in the algorithm is tournament.

MOICA
This algorithm, first proposed by Atashpaz & Lucas [24], that is the population-based algorithm, similar to the other meta-heuristic algorithms in which the space of the solution is searched by points in the name of the country (similar to the chromosome in GA). In this algorithm we deal with the imperialist and dominant countries. colonial and colonial states that make up the colonies. With colonial competition, the dominant countries are attracted to colonialism and they form an empire. Over time, the power of the countries under the control of an empire may be greater than that of its imperialist, and thus empires will change. Finally, by creating a single empire, we get the optimal solution. Also, the flowchart of the algorithm used in this paper is as follows (see Fig. 2).

Measuring Metrics
Multi-objective models need metrics to compare. Usually these criteria include four main items, which are shown below:

Parameters Setting
The performance of the meta-heuristic algorithms is usually sensitive to the settings of the parameters that can influence the search behavior. In this section, the parameters of the two algorithms are discribed. For this purpose we have used the Taguchi method. The Taguchi method is one of the most powerful statistical methods used to set parameters.

Proposed NSGA-II
In the NSGA-II, four parameters, namely population size, mutation rate, crossover rate and iteration number are used with three levels for each parameter. Given the number of factors and analysis levels, the standard orthogonal table L9 provided by Taguchi method, from MINITAB software, was chosen for this study.

Proposed MOICA
In the MOICA algorithm the following input parameters are considered: the number of countries and imperialist, assimilation coefficient, number of iterations, assimilation angle coefficient and alpha coefficient. Given the number of factors and analysis levels, the standard orthogonal table L27 provided by Taguchi method, from MINITAB software, was chosen for this study. The results are in Tab. 4.

DISCUSSION
This section discusses the results of model solving by the mentioned algorithms. For this purpose, according to the dimensions of the problem, the model is divided into three sections: small, medium and large size. Model results are also executed by a computer with cpu 2.67 GHz and 4 G RAM.

Model Results
Considering the division of sample problems, we consider three models (P 01 , P 02 and P 03 ). The uncertainty index for demand is equal to nominal value, as well as four levels of uncertainty for non-deterministic parameters (i.e., ρ = 0.2, 0.5, 0.7, 1). First, the deterministic and robust models are solved under nominal data. In Tab. 5 nominal data are generated using the random distributions (randomly). Then, under each uncertainty level (ρ), three random scenario are generated in the corresponding uncertainty set to demonstrate the performance of the solutions. In addition, the mean and standard deviation of the objective function are assesed for three scenarios at each level of uncertainty (for two algorithms). The results of experiments under customer demands are reported in Tab. 6 and Tab. 7. The results showed that the robust model obtained the solutions with both higher quality and lower standard deviations than the deterministic model. In all problems, with respect to these tables except problem P01 with ρ = 0.2 the robust method dominates the deterministic one with respect to the two relevant criteria. In addition, approximately the objective function is generally increased by increasing the level of uncertainty.
In the following, we compare the values for the objective function and the standard deviation obtained from the two algorithms with respect to the above tables. The results show that in almost all cases (sample problems P01, P 02 and P 03 ) the robust model has a better quality than the deterministic model (Except P 01 with ρ = 0.7). For example, we consider the problem P 01 , and Fig. 3 and Fig.  4 in this case show that the performance of the robust model for both mentioned meta-heuristic algorithms is better than the deterministic model. For the other problems (P 02 and P 03 ), according to the results presented in Tab. 6 and Tab. 7, the same result has been achieved.  As shown in Fig. 5, with increasing the size of the problem, the distance between the solutions of the two models also increases. This means that in larger-sized problems, the robust model offers more qualitative solution than the deterministic model. The Robust model also has the ability to solve the model with the maximum number of non-deterministic parameters according to the size of each problem. Therefore, the effectiveness of the robust model is verifiable.
Finally, we present a summary of the results of the implementation of two algorithms for the robust model. It is noteworthy that since the lower values of MID and SM metrics imply better performance of the algorithm, the alternative hypothesis of a problem is changed 0 d < . This hypothesis is tested using t-test at 97.5% level of significance. Additionally, the Kolmogorov-Smirnov test is also employed to test the collected data for normality, because t-test and ANOVA assume that the sample data follow a normal distribution. In the K-S test, it is important to note that if the p-value is greater than 0.05, the data can be normalized. From Tab. 9 we find that only with respect to the DM metric, the NSGA-II dominate the MOICA algorithm. However, the MOICA method takes less time than the NSGA-II algorithm. So in general it can be said that the NSGA-II algorithm has a better quality than the MOICA in IRP issues. To shed some light on the issue, we will, for example, consider the problem P02 (with ρ = 0.5) and compare Pareto's solutions to each other. It should be noted that from the mean values of Pareto's solutions, the values of the objective functions (Z1, Z 2 ) for each test problem are obtained. Finally, Pareto's solutions to these two algorithms are shown in Fig. 6. Also, as can be seen, the Pareto solutions of the NSGA-II are superior to the MOICA algorithm.

CONCLUSION
From the past to today, inventory has always been important to managers. For example, warehousing issues have been raised since the 60s and have been instrumental in reducing system costs. On the other hand, transportation issues are also considered for each company's goals. Moving goods at a low cost is a major factor in the cost segment. But this factor alone is not enough in companies' productivity. For this reason, in the 1990s, the combination of these two issues led to the emergence of the IRP, and from that year on, various versions of it were created. For example, IRPB can be noted that it has many applications in the real world that are also mentioned in this research.
In this paper, a new mathematical model and solution approach for a bi-objective robust inventory routing problem was proposed. The first objective function included the system costs and the second objective function minimizes transportation risks on routes taken by vehicles. Also we considered some of the parameters of the model, such as customer demand, in an uncertain form and based on robust optimization approach (bertsimas and sim, 2004) and two meta-heuristic algorithms, namely MOICA and NSGA-II, were developed to solve the proposed model. In the following, we compared the robust model with a deterministic model based on two mentioned algorithms. The results indicated that the robust model has a better quality than the deterministic model in terms of the mean Pareto's solutions of objective functions and the standard deviation of these solutions and so the efficiency of the robust approach is proven in this paper.
It should be noted that according to the literature review, the presentation of the IRPB model with multiproducts and multi-periods, in this paper, is considered to be the contribution of this research.
Concerning future research, it can be said that the development of robust models can be of interest to other authors. For example, more general forms of uncertainty sets (such as polyhedral uncertainty sets) can be considered in developing robust optimization models for inventoryrouting problems. It is also useful to use other approaches to solve non-deterministic models such as fuzzy method and its combination with other methods, as well as the development of exact methods for solving small-scale problems (e.g., branch and price or benders decomposition).