# A Modified Priority-Based Encoding for Design of a Closed-Loop Supply Chain Network Using a Discrete League Championship Algorithm

Authors: Javid Ghahremani Nahr; Ramez Kian; Hassan Rezazadeh
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In a closed-loop supply chain network, the aim is to ensure a smooth flow of materials and attaining the maximum value from returning and end-of-life goods. This paper presents a single-objective deterministic mixed integer linear programming (MILP) model for the closed-loop supply chain (CLSC) network design problem consisting of plants, collection centers, disposal centers, and customer zones. Our model minimizes the total costs comprising fixed opening cost of plants, collection, disposal centers, and transportation costs of products among the nodes. As supply chain network design problems belong to the class of NP-hard problems, a novel league championship algorithm (LCA) with a modified priority-based encoding is applied to find a near-optimal solution. New operators are defined for the LCA to search the discrete space. Numerical comparison of our proposed encoding with the existing approaches in the literature is indicative of the high quality performance of the proposed encoding.


---

## Body

## 1. Introduction

Nowadays, supply chain management (SCM) has received attention in several organizations. SCM is described as the design, production, organization, execution, control, and testing regarding supply chain activities with the goal of creating net value, minimizing the logistics cost, creating a competitive infrastructure, synchronizing demand with supply, leveraging worldwide supply chain, and measuring effectiveness globally. Generally, SCM includes the coordination and integration of key business activities including activities from purchase of raw materials to distribution of the finished products to customers. An efficient and effective supply chain can be regarded as a competitive advantage for companies and plants and helps them to cope with the global market pressure. There are two kinds of supply chain, i.e., forward supply chains and reverse supply chains. The forward supply chain is defined as a set of activities converting raw materials to products as well as storing and distributing products to the customers, while the reverse supply chain consists of a series of activities such as collection, inspection, repair, recovery, and disposal of used products. Integration of reverse and forward supply chains creates a CLSC. In other words, both forward and reverse supply chain networks are present in the CLSC networks. Network design is one of the most significant strategic decisions in SCM. In general, supply chain network design decisions include determining the number and location of facilities and the quantity of flow between them. In recent years, a few studies have focused on integrated forward and reverse network designs, while this type of integration can prevent the suboptimality and increase the level of network performance and coordination between forward and reverse processes. The present paper proposes a MILP model for design of a CLSC network consisting of plants, collection centers, disposal centers, and customer zones. Furthermore, it addresses two different problems, i.e., the facility location problem and the quantity of flow between facility optimization. A discrete LCA with a new modified priority-based encoding is applied to solve the proposed model to minimize the total network costs.


## 2. Literature Review

Over recent years, considering the rising importance of reverse and CLSC network designs, numerous articles have been published in this regard. Govindan et al. [1] present a more comprehensive literature review regarding the closed-loop and reverse supply chains. They classify 382 papers published from 2007 to 2013 and propose a more detailed classification based on 10 factors, e.g., the year of release, approaches, objectives, and functions. They assert that almost 50% of the total surveys are linked to the CLSC network design, and almost 40% of them are connected to the reverse supply chain network design (RSCND). Furthermore, the mentioned study reveals that 12% and 88% of the published papers are related to the single-objective and multiobjective models, respectively. Most of the logistics network (RL) design problems (both forward and reverse logistics) include different MILP-based facility location models. These models include vast varieties from simple ones, such as locating facilities with unlimited capacity and lightweight single-piece model, and single product (e.g., [2]) to more complex ones, such as models for limited-capacity multilevel or multiproduct and multiperiod models (e.g., [3]). Jayaraman et al. [4] develop a MILP model for RL network design under a pull system based on customer demand for recycled products aiming to minimize the total cost. Krikke et al. [5] present a MILP for a two-stage RL network of a copier production plant. In their model, both processing returns and inventory costs are considered in the objective function. The uncertainty of the return number and quality of products is an important factor in the design of the RL network. Accordingly, Listeş and Dekker [6] present a MILP model for a sand recycling network with the aim of maximizing the total profit. They extend their model to the different conditions under several scenarios. Aras et al. [7] present a nonlinear model to determine the location of collection centers in a simple RL network. The most important thing to note in their paper is the ability of the proposed model in determining optimal purchase price of the used products with the profit-maximizing objective function. In their solution approach, they develop a heuristic method based on Tabu Search (TS) algorithm. Üster et al. [8] develop a semi-integrated network in which there is a forward logistics network where only collection and recycling centers should be located. It optimizes forward and reverse flows simultaneously. They propose an exact approach based on Benders decomposition technique to solve the model. Lu and Bostel [9] consider a two-level location problem with three types of facilities that should be in a special reverse logistics called Reconstruction Location Network. They present a mixed integer binary programming model in which the forward and reverse flow and the interaction between them are considered at the same time. They also develop an algorithm based on Lagrangian heuristic algorithms to solve the model. Wojanowski et al. [10] study the interactions between industries and government agencies in relation to a series of products used by families. In order to design a CLSC for third-party logistics providers Du and Evans [3] propose an advanced biobjective MILP model. The objective functions of their model are minimizing the lateness and the total cost. They develop a hybrid scatter search method to solve their model. Pishvaee et al. [11] propose a linear model for the location of collection and inspection centers in a reverse logistics network and develop a simulated annealing (SA) algorithm to solve the model in large-scale sizes. Bing et al. [12] consider the plastic recycling in the Netherlands. They propose a MILP model and aim to minimize transportation costs and environmental impacts under different scenarios to find the optimal separation strategy. Gomes et al. [13] propose a MILP model to find the optimal location for the collection and sorting centers. Selection is carried out simultaneously under what called “Tactical Network Planning.” Their work is inspired by the European Union directive for electrical and electronic waste. Alumur et al. [14] propose a general MILP model with great flexibility that can be used for different recycled products and can be expanded further to include more settings. Their model determine optimal locations and capacities of inspection centers and reproduction plants in the design of RL network. Also, a case study with a RL network design concept for washing machines and dryers (large appliances) was studied in Germany. Their important assumption is that all the components of a product can be used again. Kannan et al. [15] propose a multistage, multiperiod, multiproduct integrated forward/reverse logistics network model for returned products, which is based on GA-based heuristic algorithm. Listeş and Dekker [6] present a MILP model and certain dynamic aspects such as due dates and inventory status which results in a complicated model and for this reason they consider a single level single product network and solved it by genetic algorithms. Giri and Sharma [16] solve the problem with developed algorithms for sequential and global optimization. Pati et al. [17] perform a different formulation by presenting a MILP model for a better management of recycled paper logistics system. The authors examine the relationship between multiple targets of recycled paper distribution network such as RL costs, improving product quality and benefits of recycling waste-paper. In addition, their model determines strategic decisions such as locating facilities, as well as tactical decisions, such as recyclable product flows and routing under multi-item, multilevel, and multifacility. Godichaud and Amodeo [18] perform a multiobjective optimization in order to increase control policies with regard to returned products. Fonseca et al. [19] consider uncertainty with the cost of transportation and the waste of production. They present a comprehensive model for reverse logistics planning that consider levels of multiple facilities, multiple products, and selection of different technologies apart from being random. They present a two-stage biobjective mixed integer stochastic programming model for providing strategic and tactical decisions, respectively, in the first and second stages to minimize costs and also the negative impact. Özceylan et al. [20] develop a mixed integer nonlinear programming (MINLP) model, which optimizes the tactical decisions on balancing the decomposition lines in the reverse supply chain and the strategic decisions related to the quantity of products flowing on the CLSC. Soleimani et al. [21] propose a stochastic integer linear programming model for designing a multiproduct CLSC in the case of uncertainties regarding the demand, purchase price, and rate of return. Some of the most crucial studies addressing the supply chain are based on a metaheuristic algorithm (Table 1). However, recently other researchers have investigated this problem meticulously (e.g., [22, 23]). Table 1 summarizes similar works and highlights our contribution in this study.Table 1
Papers published in the field of SCND using metaheuristics algorithms.
Author-Date Model Network Encoding Algorithm Gen et al. [25] MILP SC Priority-based encoding Genetic Algorithm Pishvaee et al. [11] MILP RL Priority-based encoding Simulated Annealing Costa et al. [26] MILP SC Matrix encoding Genetic Algorithm Jamshidi et al. [27] MILP SC Priority-based encoding Hybrid Memetic Algorithm Subramanian et al. [28] MILP CLSC Priority-based encoding Simulated Annealing Zohal and Soleimani [29] MILP CLSC Matrix encoding Ant Colony Algorithm Zandieh and Chensebli [30] MILP RL Priority-based encoding Water Flow-like Algorithm Singh et al. [31] MILP RL Matrix encoding Genetic Algorithm Pasandideh and Asadi [32] MILP SC Priority-based encoding Imperialist Competitive Algorithm This Paper MILP CLSC A new Priority-based encoding Discrete League Championship Algorithm


## 3. Contribution Highlights

We contribute to the literature in three dimensions: (i) to the best of our knowledge, this work is the first League Champion Algorithm in which continuous encoded solutions are converted to discrete encoded solutions based on the new operators defined in this paper. (ii) We propose an efficient encoding methodology in our algorithm which keeps the generated solutions in the feasible region which causes speeding up of the algorithm and rising the chance of finding optimal solution. (iii) The numerical experiments demonstrate that our algorithm outperform similar existing ones in the literature and capable of solving large-scale problems in reasonable execution time which are not solvable via the commercial solver in reasonable time.


## 4. Problem Definition and Modeling

In the current study, a general CLSC network is considered. The forward network includes plants and customer zones, and the reverse network includes collection and disposal centers. According to Figure1, the plants could manufacture new products and remanufacture the returned ones. The plants send products to the customer zones. Then, the returned products from customers are collected by the collection centers, and after inspection of products, the repairable products are sent to the plants and the remaining ones are sent to the disposal center. To determine the scope of the study, the following assumptions are made for the proposed model.(i)
The model is designed for a single period.(ii)
All facilities have limited and identified capacities.(iii)
The locations of all centers are potential and unknown.(iv)
All customers demands must be satisfied and all the returned products from customers must be collected.(v)
Customers locations are fixed and predefined.Figure 1
The closed-loop supply chain network Amin and Zhang [24].Given the above assumptions, the most important issue mentioned in this paper is locating the plants, collection centers, and disposal centers as well as determining the optimal amount of flow between centers. The network can be formulated as a MILP model. Sets, parameters, and decision variables are defined as follows:Sets I: set of potential plants, whose elements are addressed by index i∈I J: set of fixed locations of customer zones, whose elements are addressed by index j∈J K: set of potential collection centers, whose elements are addressed by index k∈K L: set of potential disposal centers, whose elements are addressed by index l∈LParameters T P S i j: unit transportation cost between plant i and customer zone j T S C j k: unit transportation cost between customer zone j and collection center k T C P k i: unit transportation cost between collection center k and plant i T C D k l: unit transportation cost between collection center k and disposal center l E i: fixed cost for opening plant i F k: fixed cost for opening collection center k G l: fixed cost for opening disposal center l C a p P i: capacity of plant i C a p C k: capacity of collection center k C a p D l: capacity of disposal center l D e m j: demand of customer j r j: return of customer j β: minimum disposal fractionDecision Variables X i j: quantity of new products shipped from plant i to customer zone j Y j k: quantity of returned products from customer j to collection center k S k i: quantity of returned products from collection center k to plant i T k l: quantity of returned products from collection center k to disposal center l Z i: 1, if a plant is located and set up at potential site i, and 0, otherwise W k: 1, if a collection center is located and set up at potential site k, and 0, otherwise H l: 1, if a disposal center is located and set up at potential site l, and 0, otherwise.


### 4.1. Model Formulation

The mathematical model of the problem can be presented as follows:(1)Minimize:∑iEiZi+∑kFkWk+∑lGlHl+∑i∑jTPSijXij+∑j∑kTSCjkYjk+∑k∑iTCPkiSki+∑k∑lTCDklTkl,(2)s.t.∑iXij≥Demj,∀j∈J,(3)∑iXij≥∑kYjk,∀j∈J,(4)∑kYjk=rj,∀j∈J,(5)β∑jYjk=∑lTkl,∀k∈K,(6)∑jYjk=∑iSki+∑lTkl,∀k∈K,(7)∑jXij+∑kSki≤ZiCapii,∀i∈I,(8)∑jYjk≤WkCapCk,∀k∈K,(9)∑kTkl≤HlCapDl,∀l∈L,(10)Zi,Wk,Hl∈0,1,∀i∈I,k∈K,l∈L,(11)Xij,Yjk,Ski,Tkl≥0,∀i∈I,j∈J,k∈K,l∈L.The first term in the objective function (1) represents the fixed costs of locating the plants. The second and third terms indicate the fixed costs of locating the collection centers and disposal centers, respectively. The fourth term corresponds to the production and transportation costs of new products. The fifth term represents the collection processing and transportation costs of returned products from customers. The sixth term calculates the recovery processing and transportation costs of returned products from collection centers to plants. The seventh term calculates the disposal processing and transportation costs of returned products from collection centers to disposal centers. Constraint (2) ensures that all customers demands are satisfied. Constraint (3) guarantees that forward flow is greater than reverse one. Constraint (4) computes the returned products from each customer. Constraint (5) enforces a minimum disposal fraction for each product. Constraint (6) indicates that the quantity of returned products from customer zones to collection centers is equal to the quantity of returned products from collection centers to plants and quantity of returned products from collection centers to disposal centers. Constraint (7) is a capacity constraint of plants. Constraint (8) is a capacity constraint for collection centers. Constraint (9) is a capacity constraint of disposal centers. Constraint (10) and (11) enforce the nonnegativity and binary restrictions on decision variables. The NP-hardness of supply chain network design problem has been proved by many research studies (e.g., [4]). The considered model in this paper consists of two different problems, i.e., facility location problem and quantity of flow between facility optimization; therefore, the developed model is reducible to facility location problem. Davis and Ray [33] concluded that facility location problem is NP-complete. Hence, the discussed CLSC network design problem is considered as NP-Hard in this paper. Solving this problem by exact solutions is time-consuming and sometimes impractical in large scales. Therefore, several metaheuristic algorithms with different approaches have been developed to get near optimal solutions, though all are not efficient. In the this paper, a discrete LCA solution approach is proposed and applied based on modified priority-based encoding.


## 4.1. Model Formulation

The mathematical model of the problem can be presented as follows:(1)Minimize:∑iEiZi+∑kFkWk+∑lGlHl+∑i∑jTPSijXij+∑j∑kTSCjkYjk+∑k∑iTCPkiSki+∑k∑lTCDklTkl,(2)s.t.∑iXij≥Demj,∀j∈J,(3)∑iXij≥∑kYjk,∀j∈J,(4)∑kYjk=rj,∀j∈J,(5)β∑jYjk=∑lTkl,∀k∈K,(6)∑jYjk=∑iSki+∑lTkl,∀k∈K,(7)∑jXij+∑kSki≤ZiCapii,∀i∈I,(8)∑jYjk≤WkCapCk,∀k∈K,(9)∑kTkl≤HlCapDl,∀l∈L,(10)Zi,Wk,Hl∈0,1,∀i∈I,k∈K,l∈L,(11)Xij,Yjk,Ski,Tkl≥0,∀i∈I,j∈J,k∈K,l∈L.The first term in the objective function (1) represents the fixed costs of locating the plants. The second and third terms indicate the fixed costs of locating the collection centers and disposal centers, respectively. The fourth term corresponds to the production and transportation costs of new products. The fifth term represents the collection processing and transportation costs of returned products from customers. The sixth term calculates the recovery processing and transportation costs of returned products from collection centers to plants. The seventh term calculates the disposal processing and transportation costs of returned products from collection centers to disposal centers. Constraint (2) ensures that all customers demands are satisfied. Constraint (3) guarantees that forward flow is greater than reverse one. Constraint (4) computes the returned products from each customer. Constraint (5) enforces a minimum disposal fraction for each product. Constraint (6) indicates that the quantity of returned products from customer zones to collection centers is equal to the quantity of returned products from collection centers to plants and quantity of returned products from collection centers to disposal centers. Constraint (7) is a capacity constraint of plants. Constraint (8) is a capacity constraint for collection centers. Constraint (9) is a capacity constraint of disposal centers. Constraint (10) and (11) enforce the nonnegativity and binary restrictions on decision variables. The NP-hardness of supply chain network design problem has been proved by many research studies (e.g., [4]). The considered model in this paper consists of two different problems, i.e., facility location problem and quantity of flow between facility optimization; therefore, the developed model is reducible to facility location problem. Davis and Ray [33] concluded that facility location problem is NP-complete. Hence, the discussed CLSC network design problem is considered as NP-Hard in this paper. Solving this problem by exact solutions is time-consuming and sometimes impractical in large scales. Therefore, several metaheuristic algorithms with different approaches have been developed to get near optimal solutions, though all are not efficient. In the this paper, a discrete LCA solution approach is proposed and applied based on modified priority-based encoding.


## 5. Solution Approach

In this section, we first discuss about the encoded solution format which is referred to as “representation” in the rest of this manuscript and describe the discrete LCA used for solving the CLSC network design problem will be explained.


### 5.1. Representation

Representation is one of the most essential issues for encoding and decoding, which affects the performance of algorithms. Tree-based representation is one way for representing network problems. There are several ways to encode the tree; for example, Gen and Cheng [34] introduce three ways of encoding the tree, i.e., edge-based encoding, vertex-based encoding, and edge-vertex encoding. Michalewicz et al. [35] use matrix-encoding representation, which belongs to edge-based encoding for solving the transportation tree. They present solutions through a K×J matrix in their approach, where K and J are the number of sources and depots, respectively. This method requires a large memory space on the computer. Another representation of transportation tree belonging to the vertex-based encoding is the PrÜfer number, which is developed by Gen and Cheng [34] wherein the solution is presented through |K|+|J|-2 digits. Their method needs some repair mechanisms to obtain a feasible solution. Gen and Cheng [34] develop a priority-based encoding for solving the transportation tree, which does not need any excessive repair mechanisms and the solution is presented through K+J digits. Moreover, this representation is applied to the shortest path problem and projects the scheduling problem. In the present paper, a modified priority-based encoding representation of Gen and Cheng [34] is proposed. In contrast to the other developed representations, presented in the literature, the suggested representation can solve the facility location and quantity flow optimization problems together. The representation includes four segments of sizes K+J, I+J, I+K, and K+L digits, respectively, wherein the position of each digit represents the sources and depots within the supply chain network.The first segment is devoted to the customers and collection centers, the second segment is devoted to the plants and customers, the third segment is dedicated to the collection centers and plants, and finally the fourth segment is related to the collection and disposal centers. For example, the CLSC network depicted in Figure2 includes five customers, two collection centers, three plants, and two disposal centers with the corresponding representations.Figure 2
An illustration of the CLSC network and its representation.To demonstrate the effectiveness of the proposed representation in solving the supply chain problem, Figure3 depicts the difference between the priority-based encoding developed by Gen et al. [25], priority-based encoding modified by Jamshidi et al. [27], and priority-based encoding proposed in this paper, which are denoted by (PB1), (PB2), and (PB3), respectively. This figure illustrates a two-level supply chain with 3 sources and 4 depots, their corresponding capacities, depots demand, the opening cost of sources, transportation cost between the nodes, and the priority-based encoding. In the proposed representation, the solution is presented through K+J digits, where the position of each digit represents the sources and depots within the supply chain network. Furthermore, the value in digits denotes the priorities. After the encoding operation, in order to create a connection between the representation and the supply chain network, the decoding should be performed in a specific manner in which the output denotes the facilities opened and transportation amount between opened centers, as well. To decode the representation, first the potential facilities is located; then, the optimal shipment size among the located centers is determined. For example, to decode a two-level supply chain network given in Figure 3, the following steps are taken into account:Figure 3
Samples of two-level supply chain network and its encoding.LetTrij denote the transportation cost from source i to depot j and let v(·) denote the priority.In the first step of the first section of decoding procedure:(i)
The source with the highest priority (source 2) is selected; then, the capacity of this source is compared with the total demand of depots. If the capacity of the selected source is less than total demand (950<1300), then the source with the next highest priority (source 3) is selected. Continue this procedure until the total capacity of sources is less than the total demand of depots. Then, reduce the priorities of the nodes (sources) which are not selected to zero; i.e., v(1)=0, and set the transportation cost from them to infinity, Tr1j=∞,∀j∈J.In the next step,(i)
the depot with the highest priority (depot 1) is assigned to the selected sources (source 3) as they have the lowest transportation cost among other pairs;(ii)
among the selected nodes, determine the shipment size; here, it is equal toX31=min⁡(700,300)=300;(iii)
Update the capacity and the demand of the selected source and depot asCap3=700-400=300 and Dem1=300-300=0;(iv)
As the demand of depot (1) is zero, its priority must be reduced to zero; i.e.,v(4)=0.In the next step,(i)
depot (4) with the highest priority after updating is connected to source (3), and shipment size between them is determined. At the end, capacity values of source (3) and the demand of depot (4) are updated. This sequence of operations is repeated until all demand of depots is satisfied and all priorities are reduce to zero. Table2 presents the trace table for the example two-level supply chain network given in Figure 3 to show how its modified priority-based encoding is obtained. In this table, column v(k) denotes the priority of the source nodes while v(k∗) denotes the priority values of the sources which will serve and set the others to zero. The capacity of the sources, the total capacity of the selected nodes, and demand vectors are given in columns Cap, CapT, and Dem, respectively. Column v(k∗+j) gives the representation and column K∗ nominate the selected source whereas column Xk∗j shows the flow amount from the source k∗ to the demand node j.Table 2
Trace table of decoding procedure Iteration.
Iteration v ( k ) Cap C a p T ∑j=1JDemj k v(k∗) 0 [2 5 3] ( 800,950,700 ) 950 1300 2 [0 5 0] 1 [2 0 3] ( 800,0 , 700 ) 1650 1300 3 [0 5 3]  Iteration v(k∗+j) C a p D e m K ∗ j X k ∗ j  0 [0 5 3 ∣ 7 4 1 6] ( 950,700 ) ( 300,350,300,350 ) 3 1 300 1 [0 5 3 ∣ 0 4 1 6] ( 950,400 ) ( 0,350,300,350 ) 3 4 350 2 [0 5 3 ∣ 0 4 1 0] ( 950,50 ) ( 0,350,300,0 ) 2 2 350 3 [0 0 3 ∣ 0 4 1 0] ( 600,50 ) ( 0,0 , 300,0 ) 2 2 - 4 [0 0 3 ∣ 0 0 1 0] ( 600,50 ) ( 0,0 , 300,0 ) 3 2 - 5 [0 0 0 ∣ 0 0 1 0] ( 600,50 ) ( 0,0 , 300,0 ) 2 3 300 6 [0 0 0 ∣ 0 0 0 0] ( 300,50 ) ( 0,0 , 0,0 )The decoding algorithm for representation of a two-level supply chain network is indicated in Algorithm1. The representation should be divided to four segments to apply the above-presented decoding algorithm to the discussed CLSC network design.Algorithm 1:The decoding algorithm of each section of representation for two-level supply chain network.
Require: Gets K~: Set of sources, J~: Set of depots, Demj: demand on depot j, Capk: capacity of source k Ensure: gives U~kj: Quantity of shipment between source k and depot j (1)CapT:=0, K′=K, K′′=∅ (2)while  CapT<∑j=1JDemj  do (3)select a node on n=argmax⁡{v(k),k∈K1′} (4)K′=K′∖{n} (5)K′′=K′′∪{n} (6)CapT=∑n∈K′Capn (7)end while (8)v(n)=0,∀n∈K′ (9)Trnj=∞,∀j∈J~,n∈K′ (10)Xkj=0,∀j∈J,k∈K (11)while  v(n)≠0,∀n∈J∪K~  do (12)Select a node based on n=argmax⁡{v(t),t∈K~+J~} (13)if  l∈k a source is selected k∗=n  then (14)j∗=argmin⁡{Trkj},j∈J~ select a depot with minimum transportation cost (15)elseIf  n∈J a depot is selected j∗=n (16)k∗=argmin⁡{Trkj∣k∈K′′} select a source with minimum transportation cost (17)end if (18)U~k∗j∗=min⁡(Capk∗,Demj∗) (19)Update demand and capacities: (20)Capk∗=Capk∗-U~k∗j∗ (21)Demj∗=Demj∗-U~k∗j∗ (22)If  Capk∗=0  then  v(k∗)=0 (23)If  Demj∗=0  then  v(j∗)=0 (24)end whileDecoding of the second segment of the representation is impossible until the first segment is decoded because customers demands consist of new products and recovered products, and the amount of recovered products is achieved only by decoding the first segment. Then, decoding of the third and fourth segments starts after determining the number and location of collection centers and the amount of the returned products. Decoding algorithm of representation for multilevel CLSC network is provided in Algorithm2.Algorithm 2:The decoding algorithm of representation for multilevel CLSC network.
Requires: problem parameters Ensure: Calculating the objective function in the multi-level CLSC (1)Zi=0,Wk=0;Hl=0,Xij=0,Yjk=0,Ski=0,Tkl=0,∀i∈I,j∈J,k∈K,l∈L (2)Calculate  Yjk≔U~jk,∀j∈J,k∈K  by Algorithm 1 (3)If  ∑jYjk>0  then  Wk=1 (4)Calculate  Tkl≔U~kl,∀k∈K,l∈L  by Algorithm 1 (5)If  ∑kTkl>0  then  Hl=1 (6)Calculate  Xij≔U~ij,∀i∈I,j∈J  by Algorithm 1 (7)If  ∑jXij>0  then  Zi=1 (8)Calculate  Ski,∀k∈K,i∈I  by Algorithm 1 (9)Calculate the total cost: (10)obj=∑i=1IEiZi+∑k=1KfkWk+∑l=aLGlHl+∑i=1I∑j=1JTPSijXij (11)+∑j=1J∑k=1KTSCjkXjk+∑k=1K∑i=1ITCPkiXki+∑k=1K∑l=1LTCDklXkl


### 5.2. Discrete League Championship Algorithm

The LCA, first presented by Kashan [36], is a population-based algorithm for global search in a continuous space, which is inspired by the championship process of sports leagues in the real world. One of the common characteristics shared by all population-based algorithms such as LCA is their attempt to move a population of possible solutions to a number of promising areas in search of a desirable solution. Similar to many population-based algorithms, a set of l solutions from the search space is randomly selected, which constitutes the initial population in the algorithm. Each solution obtained from the population is related to one of l teams (l is an even number) that shows the current formation of the team. Hence, team i represents ith member of the population. Each solution obtained from the population has its own fitness value. In this algorithm, different solutions can be provided for a problem, which are compared on the basis of their fitness values (their objective function values). Each of the solutions is improved and ultimately near-optimal solution is selected. A number of teams (examined random initial solutions) compete together in a pairwise form as a league within a few weeks (the number of evaluation steps in an iterative algorithm). Winner and loser teams are determined (draw is not allowed) based on their power of play (metaphor of the fitness or objective value of the solution) resulting from the team formation. Every week, each team forms a new team formation by its coach with an artificial process analyzing the last week’s matches and obtains the best recognized formation up to that time (iterations number of the algorithm), and this process will continue. A schematic overview of our algorithm is given in Figure 4.Figure 4
The overview of the LCA algorithm implementation.ParametersS (examined number of seasons) and l (the number of teams) and constant coefficients used to scale the contribution of the strength and weakness of components (denoted by ψ1 and ψ2), are adjustable parameters whose changes have a direct impact on the final answer of the algorithm. The solution space in the CLSC problem is discrete, which means that components of each individual solution within the population cannot get an arbitrary amount, and they are allowed only to get natural values. Hence, a discrete LCA is dealt with to solve the CLSC problem. In addition, the numbers in representations should not be repeated. In the conducted studies dealing with discrete league championship algorithm, the algorithm searches a continuous space and finally presents a discrete solution by an innovative method. However, in the procedure proposed in the present study, all stages of the algorithm are made in the discrete space with no changes in the presentation of solution.Hence, according to the discrete structure of the problem and the definition of team formation, discrete version of the operators (including addition, multiplication, subtraction, and arranging) is defined to address the strategies used in LCA in an appropriate structure to solve the CLSC problem. The mentioned operators and their corresponding arguments are summarized below:(1)
Arrange (formation, strategy)→Plus Formation(2)
Subtraction (formation, formation)→Minus Strategy(3)
Addition (strategy, strategy)→Plus Strategy(4)
Multiplication (real_number, strategy)→Times StrategyThe strategy vector and the above mentioned operators will be mathematically defined.


#### 5.2.1. Strategy Vector and Organization Operator (Arrange)

The strategy vector, denoted byS, should be defined in such a way that by applying to a team formation vector at any time step, another team formation is achieved. Thus, the strategy vector S for each team is defined as a sequence of “displacements” of team formation components as follows: (12)S=i1,j1,i2,j2,…,ik,jk,…,iS,jS, where S specifies the length of sequence. Applying strategy S to a vector of team formation means that the components of i1 and j1 in the vector of team formation are “displaced” together and then their components i2 and j2 and i3 and j3 are “displaced” together up to iS and jS. Hence, organization operator, i.e., the sum of team formation with a strategy, is defined as a set of displacements specified by the strategy vector in the vector of team formation. For example, if F=(1,2,3,4,5) is the vector of team formation and S=((1,2),(2,3)) is the vector of strategy, then the following situations are achieved, respectively, by applying S to F:Displacement (1) with (2):(2,1,3,4,5)Displacement (2) with (3):(3,1,2,4,5)Hence, arrange(F,S)=F+S=(3,1,2,4,5)Symmetry of a strategy vectorS is displayed as -S, which means that the sequence of displacements in S is considered reversely.


#### 5.2.2. Subtraction Operator

Ifx1 and x2 are two vectors of a team formation, the difference between x1 and x2 is defined in a way that x2-x1 is equal to the vector of strategy, S. If it is applied to x1, x2 is achieved. The algorithm which calculates the difference between x1 and x2 should be chosen with consideration of(13)x2-x1=-x2-x1x1=x2⟶x2-x1=∅Value of zero is displayed with ∅ and is equal to a null sequence.


#### 5.2.3. Addition Operator

IfS1 and S2 are two strategy vectors, sum of them is displayed as S1⊕S2 and equals to the strategy vector which is achieved by connecting the displacement sequence of S2 to the end of sequence in S1. As some displacements cancel each others effects, the obtained sequence from connecting S2 to S1 can be smaller; for example, if S1=((1,2),(2,5)) and S2=(2,5), then we have S1⊕S2=((1,2),(2,5),(2,5)).However, since two consecutive displacements(2,5) cancel each others effects, they can be deleted; hence, S1⊕S2=(1,2). Therefore, in general, we have S1⊕S2≤S1+S2. In addition, for each strategy S, we will have S⊕-S=∅.


#### 5.2.4. Multiplication Operator

Ifc is a real number and S is a strategy vector, their multiplication, cS, is defined as follows depending on the value of c.(1) Case of0≤c<1: in this case, it is assumed that cS is equal to the integral part of the number cS; vector S is truncated so that its length is equal to cS: (14)cS=ik,jk,k=1,2,…,cS. In the special case of c=0, we have cs=∅.(2) Case ofc≥1: in this case, we have c=k+c′, where k=c is a natural number and c′=c-c∈[0,1]. Therefore, cS is defined as follows: (15)cS=S⊕S⋯S︸ktimes⊕c′S.(3) Case ofc<0: in this case, the multiplication is converted to one of the above-presented modes by writing cS as (-c)(-S).We have four strategies whose notations are inspired from SWOT (Strength, Weakness, Opportunity, Threats) analysis in strategic planning literature.S/O corresponds to the aggressive strategy in which the team is strong and the opponent is week while W/T, in contrast, corresponds to a defensive strategy in which the team is in a weak mode and the opponent is strong (threat). Two other intermediate strategies are denoted by S/T and W/O corresponding toboth strong andboth weak modes, respectively. Let Xit+1 denote the formation of team i in week t and Bit be the best formation for team i in week t. The constant coefficients ψ1 and ψ2 are used to scale the strength and weakness of components, respectively. r1 and r2 are uniform random numbers from the interval [0,1].With the application of the defined operators, equations related to strategies in discrete LCA algorithm will be written as follows.(i)
If both teamsi and l in week t win their games against teams j and k, then formation of team i in week t+1 with strategy S/T will be as follows: (16)Xit+1=Bit+ψ1r1Bit-Bkt⊕ψ1r2Bit-Bjt.(ii)
If teami wins team j and team l loses against team k in week t, team formation of i in week t+1 with strategy S/O will be as follows: (17)Xit+1=Bit+ψ2r1Bkt-Bit⊕ψ1r2Bit-Bjt.(iii)
If both teamsi and l lose their games against teams j and k in week t, team formation of i in week t+1 with strategy W/T will be as follows: (18)Xit+1=Bit+ψ1r1Bkt-Bit⊕ψ2r2Bit-Bjt.(iv)
If both teamsi and l lose their games against teams j and k in week t, team formation i in week t+1 with strategy W/O will be as follows: (19)Xit+1=Bit+ψ2r1Bkt-Bit⊕ψ2r2Bjt-Bit.


## 5.1. Representation

Representation is one of the most essential issues for encoding and decoding, which affects the performance of algorithms. Tree-based representation is one way for representing network problems. There are several ways to encode the tree; for example, Gen and Cheng [34] introduce three ways of encoding the tree, i.e., edge-based encoding, vertex-based encoding, and edge-vertex encoding. Michalewicz et al. [35] use matrix-encoding representation, which belongs to edge-based encoding for solving the transportation tree. They present solutions through a K×J matrix in their approach, where K and J are the number of sources and depots, respectively. This method requires a large memory space on the computer. Another representation of transportation tree belonging to the vertex-based encoding is the PrÜfer number, which is developed by Gen and Cheng [34] wherein the solution is presented through |K|+|J|-2 digits. Their method needs some repair mechanisms to obtain a feasible solution. Gen and Cheng [34] develop a priority-based encoding for solving the transportation tree, which does not need any excessive repair mechanisms and the solution is presented through K+J digits. Moreover, this representation is applied to the shortest path problem and projects the scheduling problem. In the present paper, a modified priority-based encoding representation of Gen and Cheng [34] is proposed. In contrast to the other developed representations, presented in the literature, the suggested representation can solve the facility location and quantity flow optimization problems together. The representation includes four segments of sizes K+J, I+J, I+K, and K+L digits, respectively, wherein the position of each digit represents the sources and depots within the supply chain network.The first segment is devoted to the customers and collection centers, the second segment is devoted to the plants and customers, the third segment is dedicated to the collection centers and plants, and finally the fourth segment is related to the collection and disposal centers. For example, the CLSC network depicted in Figure2 includes five customers, two collection centers, three plants, and two disposal centers with the corresponding representations.Figure 2
An illustration of the CLSC network and its representation.To demonstrate the effectiveness of the proposed representation in solving the supply chain problem, Figure3 depicts the difference between the priority-based encoding developed by Gen et al. [25], priority-based encoding modified by Jamshidi et al. [27], and priority-based encoding proposed in this paper, which are denoted by (PB1), (PB2), and (PB3), respectively. This figure illustrates a two-level supply chain with 3 sources and 4 depots, their corresponding capacities, depots demand, the opening cost of sources, transportation cost between the nodes, and the priority-based encoding. In the proposed representation, the solution is presented through K+J digits, where the position of each digit represents the sources and depots within the supply chain network. Furthermore, the value in digits denotes the priorities. After the encoding operation, in order to create a connection between the representation and the supply chain network, the decoding should be performed in a specific manner in which the output denotes the facilities opened and transportation amount between opened centers, as well. To decode the representation, first the potential facilities is located; then, the optimal shipment size among the located centers is determined. For example, to decode a two-level supply chain network given in Figure 3, the following steps are taken into account:Figure 3
Samples of two-level supply chain network and its encoding.LetTrij denote the transportation cost from source i to depot j and let v(·) denote the priority.In the first step of the first section of decoding procedure:(i)
The source with the highest priority (source 2) is selected; then, the capacity of this source is compared with the total demand of depots. If the capacity of the selected source is less than total demand (950<1300), then the source with the next highest priority (source 3) is selected. Continue this procedure until the total capacity of sources is less than the total demand of depots. Then, reduce the priorities of the nodes (sources) which are not selected to zero; i.e., v(1)=0, and set the transportation cost from them to infinity, Tr1j=∞,∀j∈J.In the next step,(i)
the depot with the highest priority (depot 1) is assigned to the selected sources (source 3) as they have the lowest transportation cost among other pairs;(ii)
among the selected nodes, determine the shipment size; here, it is equal toX31=min⁡(700,300)=300;(iii)
Update the capacity and the demand of the selected source and depot asCap3=700-400=300 and Dem1=300-300=0;(iv)
As the demand of depot (1) is zero, its priority must be reduced to zero; i.e.,v(4)=0.In the next step,(i)
depot (4) with the highest priority after updating is connected to source (3), and shipment size between them is determined. At the end, capacity values of source (3) and the demand of depot (4) are updated. This sequence of operations is repeated until all demand of depots is satisfied and all priorities are reduce to zero. Table2 presents the trace table for the example two-level supply chain network given in Figure 3 to show how its modified priority-based encoding is obtained. In this table, column v(k) denotes the priority of the source nodes while v(k∗) denotes the priority values of the sources which will serve and set the others to zero. The capacity of the sources, the total capacity of the selected nodes, and demand vectors are given in columns Cap, CapT, and Dem, respectively. Column v(k∗+j) gives the representation and column K∗ nominate the selected source whereas column Xk∗j shows the flow amount from the source k∗ to the demand node j.Table 2
Trace table of decoding procedure Iteration.
Iteration v ( k ) Cap C a p T ∑j=1JDemj k v(k∗) 0 [2 5 3] ( 800,950,700 ) 950 1300 2 [0 5 0] 1 [2 0 3] ( 800,0 , 700 ) 1650 1300 3 [0 5 3]  Iteration v(k∗+j) C a p D e m K ∗ j X k ∗ j  0 [0 5 3 ∣ 7 4 1 6] ( 950,700 ) ( 300,350,300,350 ) 3 1 300 1 [0 5 3 ∣ 0 4 1 6] ( 950,400 ) ( 0,350,300,350 ) 3 4 350 2 [0 5 3 ∣ 0 4 1 0] ( 950,50 ) ( 0,350,300,0 ) 2 2 350 3 [0 0 3 ∣ 0 4 1 0] ( 600,50 ) ( 0,0 , 300,0 ) 2 2 - 4 [0 0 3 ∣ 0 0 1 0] ( 600,50 ) ( 0,0 , 300,0 ) 3 2 - 5 [0 0 0 ∣ 0 0 1 0] ( 600,50 ) ( 0,0 , 300,0 ) 2 3 300 6 [0 0 0 ∣ 0 0 0 0] ( 300,50 ) ( 0,0 , 0,0 )The decoding algorithm for representation of a two-level supply chain network is indicated in Algorithm1. The representation should be divided to four segments to apply the above-presented decoding algorithm to the discussed CLSC network design.Algorithm 1:The decoding algorithm of each section of representation for two-level supply chain network.
Require: Gets K~: Set of sources, J~: Set of depots, Demj: demand on depot j, Capk: capacity of source k Ensure: gives U~kj: Quantity of shipment between source k and depot j (1)CapT:=0, K′=K, K′′=∅ (2)while  CapT<∑j=1JDemj  do (3)select a node on n=argmax⁡{v(k),k∈K1′} (4)K′=K′∖{n} (5)K′′=K′′∪{n} (6)CapT=∑n∈K′Capn (7)end while (8)v(n)=0,∀n∈K′ (9)Trnj=∞,∀j∈J~,n∈K′ (10)Xkj=0,∀j∈J,k∈K (11)while  v(n)≠0,∀n∈J∪K~  do (12)Select a node based on n=argmax⁡{v(t),t∈K~+J~} (13)if  l∈k a source is selected k∗=n  then (14)j∗=argmin⁡{Trkj},j∈J~ select a depot with minimum transportation cost (15)elseIf  n∈J a depot is selected j∗=n (16)k∗=argmin⁡{Trkj∣k∈K′′} select a source with minimum transportation cost (17)end if (18)U~k∗j∗=min⁡(Capk∗,Demj∗) (19)Update demand and capacities: (20)Capk∗=Capk∗-U~k∗j∗ (21)Demj∗=Demj∗-U~k∗j∗ (22)If  Capk∗=0  then  v(k∗)=0 (23)If  Demj∗=0  then  v(j∗)=0 (24)end whileDecoding of the second segment of the representation is impossible until the first segment is decoded because customers demands consist of new products and recovered products, and the amount of recovered products is achieved only by decoding the first segment. Then, decoding of the third and fourth segments starts after determining the number and location of collection centers and the amount of the returned products. Decoding algorithm of representation for multilevel CLSC network is provided in Algorithm2.Algorithm 2:The decoding algorithm of representation for multilevel CLSC network.
Requires: problem parameters Ensure: Calculating the objective function in the multi-level CLSC (1)Zi=0,Wk=0;Hl=0,Xij=0,Yjk=0,Ski=0,Tkl=0,∀i∈I,j∈J,k∈K,l∈L (2)Calculate  Yjk≔U~jk,∀j∈J,k∈K  by Algorithm 1 (3)If  ∑jYjk>0  then  Wk=1 (4)Calculate  Tkl≔U~kl,∀k∈K,l∈L  by Algorithm 1 (5)If  ∑kTkl>0  then  Hl=1 (6)Calculate  Xij≔U~ij,∀i∈I,j∈J  by Algorithm 1 (7)If  ∑jXij>0  then  Zi=1 (8)Calculate  Ski,∀k∈K,i∈I  by Algorithm 1 (9)Calculate the total cost: (10)obj=∑i=1IEiZi+∑k=1KfkWk+∑l=aLGlHl+∑i=1I∑j=1JTPSijXij (11)+∑j=1J∑k=1KTSCjkXjk+∑k=1K∑i=1ITCPkiXki+∑k=1K∑l=1LTCDklXkl


## 5.2. Discrete League Championship Algorithm

The LCA, first presented by Kashan [36], is a population-based algorithm for global search in a continuous space, which is inspired by the championship process of sports leagues in the real world. One of the common characteristics shared by all population-based algorithms such as LCA is their attempt to move a population of possible solutions to a number of promising areas in search of a desirable solution. Similar to many population-based algorithms, a set of l solutions from the search space is randomly selected, which constitutes the initial population in the algorithm. Each solution obtained from the population is related to one of l teams (l is an even number) that shows the current formation of the team. Hence, team i represents ith member of the population. Each solution obtained from the population has its own fitness value. In this algorithm, different solutions can be provided for a problem, which are compared on the basis of their fitness values (their objective function values). Each of the solutions is improved and ultimately near-optimal solution is selected. A number of teams (examined random initial solutions) compete together in a pairwise form as a league within a few weeks (the number of evaluation steps in an iterative algorithm). Winner and loser teams are determined (draw is not allowed) based on their power of play (metaphor of the fitness or objective value of the solution) resulting from the team formation. Every week, each team forms a new team formation by its coach with an artificial process analyzing the last week’s matches and obtains the best recognized formation up to that time (iterations number of the algorithm), and this process will continue. A schematic overview of our algorithm is given in Figure 4.Figure 4
The overview of the LCA algorithm implementation.ParametersS (examined number of seasons) and l (the number of teams) and constant coefficients used to scale the contribution of the strength and weakness of components (denoted by ψ1 and ψ2), are adjustable parameters whose changes have a direct impact on the final answer of the algorithm. The solution space in the CLSC problem is discrete, which means that components of each individual solution within the population cannot get an arbitrary amount, and they are allowed only to get natural values. Hence, a discrete LCA is dealt with to solve the CLSC problem. In addition, the numbers in representations should not be repeated. In the conducted studies dealing with discrete league championship algorithm, the algorithm searches a continuous space and finally presents a discrete solution by an innovative method. However, in the procedure proposed in the present study, all stages of the algorithm are made in the discrete space with no changes in the presentation of solution.Hence, according to the discrete structure of the problem and the definition of team formation, discrete version of the operators (including addition, multiplication, subtraction, and arranging) is defined to address the strategies used in LCA in an appropriate structure to solve the CLSC problem. The mentioned operators and their corresponding arguments are summarized below:(1)
Arrange (formation, strategy)→Plus Formation(2)
Subtraction (formation, formation)→Minus Strategy(3)
Addition (strategy, strategy)→Plus Strategy(4)
Multiplication (real_number, strategy)→Times StrategyThe strategy vector and the above mentioned operators will be mathematically defined.


### 5.2.1. Strategy Vector and Organization Operator (Arrange)

The strategy vector, denoted byS, should be defined in such a way that by applying to a team formation vector at any time step, another team formation is achieved. Thus, the strategy vector S for each team is defined as a sequence of “displacements” of team formation components as follows: (12)S=i1,j1,i2,j2,…,ik,jk,…,iS,jS, where S specifies the length of sequence. Applying strategy S to a vector of team formation means that the components of i1 and j1 in the vector of team formation are “displaced” together and then their components i2 and j2 and i3 and j3 are “displaced” together up to iS and jS. Hence, organization operator, i.e., the sum of team formation with a strategy, is defined as a set of displacements specified by the strategy vector in the vector of team formation. For example, if F=(1,2,3,4,5) is the vector of team formation and S=((1,2),(2,3)) is the vector of strategy, then the following situations are achieved, respectively, by applying S to F:Displacement (1) with (2):(2,1,3,4,5)Displacement (2) with (3):(3,1,2,4,5)Hence, arrange(F,S)=F+S=(3,1,2,4,5)Symmetry of a strategy vectorS is displayed as -S, which means that the sequence of displacements in S is considered reversely.


### 5.2.2. Subtraction Operator

Ifx1 and x2 are two vectors of a team formation, the difference between x1 and x2 is defined in a way that x2-x1 is equal to the vector of strategy, S. If it is applied to x1, x2 is achieved. The algorithm which calculates the difference between x1 and x2 should be chosen with consideration of(13)x2-x1=-x2-x1x1=x2⟶x2-x1=∅Value of zero is displayed with ∅ and is equal to a null sequence.


### 5.2.3. Addition Operator

IfS1 and S2 are two strategy vectors, sum of them is displayed as S1⊕S2 and equals to the strategy vector which is achieved by connecting the displacement sequence of S2 to the end of sequence in S1. As some displacements cancel each others effects, the obtained sequence from connecting S2 to S1 can be smaller; for example, if S1=((1,2),(2,5)) and S2=(2,5), then we have S1⊕S2=((1,2),(2,5),(2,5)).However, since two consecutive displacements(2,5) cancel each others effects, they can be deleted; hence, S1⊕S2=(1,2). Therefore, in general, we have S1⊕S2≤S1+S2. In addition, for each strategy S, we will have S⊕-S=∅.


### 5.2.4. Multiplication Operator

Ifc is a real number and S is a strategy vector, their multiplication, cS, is defined as follows depending on the value of c.(1) Case of0≤c<1: in this case, it is assumed that cS is equal to the integral part of the number cS; vector S is truncated so that its length is equal to cS: (14)cS=ik,jk,k=1,2,…,cS. In the special case of c=0, we have cs=∅.(2) Case ofc≥1: in this case, we have c=k+c′, where k=c is a natural number and c′=c-c∈[0,1]. Therefore, cS is defined as follows: (15)cS=S⊕S⋯S︸ktimes⊕c′S.(3) Case ofc<0: in this case, the multiplication is converted to one of the above-presented modes by writing cS as (-c)(-S).We have four strategies whose notations are inspired from SWOT (Strength, Weakness, Opportunity, Threats) analysis in strategic planning literature.S/O corresponds to the aggressive strategy in which the team is strong and the opponent is week while W/T, in contrast, corresponds to a defensive strategy in which the team is in a weak mode and the opponent is strong (threat). Two other intermediate strategies are denoted by S/T and W/O corresponding toboth strong andboth weak modes, respectively. Let Xit+1 denote the formation of team i in week t and Bit be the best formation for team i in week t. The constant coefficients ψ1 and ψ2 are used to scale the strength and weakness of components, respectively. r1 and r2 are uniform random numbers from the interval [0,1].With the application of the defined operators, equations related to strategies in discrete LCA algorithm will be written as follows.(i)
If both teamsi and l in week t win their games against teams j and k, then formation of team i in week t+1 with strategy S/T will be as follows: (16)Xit+1=Bit+ψ1r1Bit-Bkt⊕ψ1r2Bit-Bjt.(ii)
If teami wins team j and team l loses against team k in week t, team formation of i in week t+1 with strategy S/O will be as follows: (17)Xit+1=Bit+ψ2r1Bkt-Bit⊕ψ1r2Bit-Bjt.(iii)
If both teamsi and l lose their games against teams j and k in week t, team formation of i in week t+1 with strategy W/T will be as follows: (18)Xit+1=Bit+ψ1r1Bkt-Bit⊕ψ2r2Bit-Bjt.(iv)
If both teamsi and l lose their games against teams j and k in week t, team formation i in week t+1 with strategy W/O will be as follows: (19)Xit+1=Bit+ψ2r1Bkt-Bit⊕ψ2r2Bjt-Bit.


## 5.2.1. Strategy Vector and Organization Operator (Arrange)

The strategy vector, denoted byS, should be defined in such a way that by applying to a team formation vector at any time step, another team formation is achieved. Thus, the strategy vector S for each team is defined as a sequence of “displacements” of team formation components as follows: (12)S=i1,j1,i2,j2,…,ik,jk,…,iS,jS, where S specifies the length of sequence. Applying strategy S to a vector of team formation means that the components of i1 and j1 in the vector of team formation are “displaced” together and then their components i2 and j2 and i3 and j3 are “displaced” together up to iS and jS. Hence, organization operator, i.e., the sum of team formation with a strategy, is defined as a set of displacements specified by the strategy vector in the vector of team formation. For example, if F=(1,2,3,4,5) is the vector of team formation and S=((1,2),(2,3)) is the vector of strategy, then the following situations are achieved, respectively, by applying S to F:Displacement (1) with (2):(2,1,3,4,5)Displacement (2) with (3):(3,1,2,4,5)Hence, arrange(F,S)=F+S=(3,1,2,4,5)Symmetry of a strategy vectorS is displayed as -S, which means that the sequence of displacements in S is considered reversely.


## 5.2.2. Subtraction Operator

Ifx1 and x2 are two vectors of a team formation, the difference between x1 and x2 is defined in a way that x2-x1 is equal to the vector of strategy, S. If it is applied to x1, x2 is achieved. The algorithm which calculates the difference between x1 and x2 should be chosen with consideration of(13)x2-x1=-x2-x1x1=x2⟶x2-x1=∅Value of zero is displayed with ∅ and is equal to a null sequence.


## 5.2.3. Addition Operator

IfS1 and S2 are two strategy vectors, sum of them is displayed as S1⊕S2 and equals to the strategy vector which is achieved by connecting the displacement sequence of S2 to the end of sequence in S1. As some displacements cancel each others effects, the obtained sequence from connecting S2 to S1 can be smaller; for example, if S1=((1,2),(2,5)) and S2=(2,5), then we have S1⊕S2=((1,2),(2,5),(2,5)).However, since two consecutive displacements(2,5) cancel each others effects, they can be deleted; hence, S1⊕S2=(1,2). Therefore, in general, we have S1⊕S2≤S1+S2. In addition, for each strategy S, we will have S⊕-S=∅.


## 5.2.4. Multiplication Operator

Ifc is a real number and S is a strategy vector, their multiplication, cS, is defined as follows depending on the value of c.(1) Case of0≤c<1: in this case, it is assumed that cS is equal to the integral part of the number cS; vector S is truncated so that its length is equal to cS: (14)cS=ik,jk,k=1,2,…,cS. In the special case of c=0, we have cs=∅.(2) Case ofc≥1: in this case, we have c=k+c′, where k=c is a natural number and c′=c-c∈[0,1]. Therefore, cS is defined as follows: (15)cS=S⊕S⋯S︸ktimes⊕c′S.(3) Case ofc<0: in this case, the multiplication is converted to one of the above-presented modes by writing cS as (-c)(-S).We have four strategies whose notations are inspired from SWOT (Strength, Weakness, Opportunity, Threats) analysis in strategic planning literature.S/O corresponds to the aggressive strategy in which the team is strong and the opponent is week while W/T, in contrast, corresponds to a defensive strategy in which the team is in a weak mode and the opponent is strong (threat). Two other intermediate strategies are denoted by S/T and W/O corresponding toboth strong andboth weak modes, respectively. Let Xit+1 denote the formation of team i in week t and Bit be the best formation for team i in week t. The constant coefficients ψ1 and ψ2 are used to scale the strength and weakness of components, respectively. r1 and r2 are uniform random numbers from the interval [0,1].With the application of the defined operators, equations related to strategies in discrete LCA algorithm will be written as follows.(i)
If both teamsi and l in week t win their games against teams j and k, then formation of team i in week t+1 with strategy S/T will be as follows: (16)Xit+1=Bit+ψ1r1Bit-Bkt⊕ψ1r2Bit-Bjt.(ii)
If teami wins team j and team l loses against team k in week t, team formation of i in week t+1 with strategy S/O will be as follows: (17)Xit+1=Bit+ψ2r1Bkt-Bit⊕ψ1r2Bit-Bjt.(iii)
If both teamsi and l lose their games against teams j and k in week t, team formation of i in week t+1 with strategy W/T will be as follows: (18)Xit+1=Bit+ψ1r1Bkt-Bit⊕ψ2r2Bit-Bjt.(iv)
If both teamsi and l lose their games against teams j and k in week t, team formation i in week t+1 with strategy W/O will be as follows: (19)Xit+1=Bit+ψ2r1Bkt-Bit⊕ψ2r2Bjt-Bit.


## 6. Numerical Results

### 6.1. Problem Instances

In order to assess the performance of the proposed discrete LCA in terms of the objective-function value and CPU time, several numerical experiments with different problem sizes were implemented, and the obtained results are reported in this section. To this end, 21 sample problems with 3 different levels, which are small, medium, and large-scale problems, were generated by different combinations of the parameter values. The levels of the generated sample problems are shown in Table3, and the ranges of the parameters are presented in Table 4. Furthermore, all parameters of the sample problems were randomly generated based on uniform distributions in prespecified intervals.Table 3
Size and scale of sample problems.
Scale of the problem Problem no. # the potential locations for plants # the potential locations for collection centers # the potential locations for disposal centers # the potential locations for customer zones Small 1 3 2 2 10 2 4 3 3 15 3 5 4 3 20 4 7 5 4 25 5 8 6 5 30 6 10 8 7 35 7 12 8 7 40  Medium 8 15 9 8 45 9 18 10 9 50 10 21 12 10 55 11 25 15 12 60 12 30 20 15 70 13 35 25 18 80 14 40 32 20 90  Large 15 45 37 24 100 16 50 42 29 120 17 60 48 33 140 18 70 53 38 160 19 80 57 42 180 20 90 60 45 200 21 100 65 50 220Table 4
Pre-specified intervals to generate parameters intervals based on uniform distributions.
Parameter Range Parameter Range Ei U ( 1000000,1200000 ) TSCjk U(20,30) Fk U ( 1000000,1200000 ) TCPki U(20,30) Gl U ( 1000000,1200000 ) TCDkl U(20,30) D e m j U(100,150) C a p P i U(800,1200) rj U ( 10,50 ) CapCk U(200,400) β U(0.6,0.8) CapDl U(100,300) T P S i j U(20,30)As the acquired results obtained from the proposed algorithm are sensitive to their initial parameters, any small changes can affect the accuracy of the best obtained solution. Therefore, the Taguchi tuning method is used for the parameters to find the best solutions. In this method, first, the appropriate factors (initial parameters) are determined, the level of each factor is selected, and then the design of experiments for this control factor is specified. After specifying the experimental design, the proposed algorithm is used to find the best combination of factors. In our study, 4 levels are considered for each factor. In Table5, the number of teams (l), number of seasons (S), and scaling constants of strength and weakness components (ψ1 and ψ2) are given.Table 5
Proposed parameter levels.
Parameter Level 1 Level 2 Level 3 Level 4 l 100 150 200 25 S 20 26 32 40 ψ1 0 2 4 6 ψ2 0 2 4 6


### 6.2. Objective Function Average

With respect to the proposed algorithm, the experimental design is employed according to the number of factors and number of levels. To this aim, each problem instance was solved five times to form the replications. Average results are reported as the final value. The best values of the proposed parameters for discrete LCA according to the mean normalized objective values are 150, 20, 6, and 2 for the number of teams, number of seasons, and constant coefficients to scale the contribution of the strength and weakness components, respectively. Moreover, Figure5 presents the normalized average of means and the S/N (Signal-to-noise) ratio graph for the experimental design of the algorithm. Based on the average of the means given in the graph, the algorithms result in a more efficient response if the parameter is placed in the lower level. Furthermore, in the S/N ratio graph, the algorithms result in a more efficient response if the parameter is placed in the higher level. According to Table 5 and Figure 5 the best performance of our algorithm is observed when the number of teams (l) and the weakness scaling coefficient (ψ2) are in level 2; the number of seasons (S) is placed at level 1, and the strength scaling coefficient (ψ4) is placed at level 4.Figure 5
Means of means and the S/N ratio plot for discrete LCA.


### 6.3. Solutions Structure and Quality

After tuning the parameters of the proposed algorithm, PB3 results are compared with those of the CPLEX outputs which were run on GAMS software (exact solution), PB2, and PB1 methods. For this purpose, five problem instances with random data given in Table4 for each problem size given in Table 3 were generated, and they were solved by each method. The average of five problems were selected as the benchmark for each sample problem. In Table 6, the average results of each sample problem instance is provided and compared with the solver’s (exact) solution.Table 6
The average of the objective values for each sample problem instance.
Scale of the problem Problem instance CPLEX PB1 PB2 PB3 Small 1 5216131.55 5216258.95 5216878.01 5216131.55 2 6537367.23 6537605.34 6538131.72 6537367.23 3 7709177.89 7709634.54 7710523.47 7709177.89 4 9378103.61 9378752.72 9379025.74 9378103.61 5 11611170.7 11612119.3 11612637.7 11611172.3 6 12217294.9 12219232.9 12218828.6 12217296.1 7 14546224.7 14547451.7 14547694.5 14546227.6  Medium 8 15642700.4 15645397.6 15644410.2 15642705.3 9 16394986.2 16411886.8 16408650.5 16395010.8 10 18192348.5 18209196.2 18199012.0 18192397.3 11 19613086.9 19800951.6 19628863.7 19613153.8 12 22984000.0 23693969.4 23485221.2 22984540.7 13 25164851.6 26212668.5 26139382.5 25165442.0 14 32156898.7 34013530.3 33511564.3 32157659.6According to the results shown in Table6, it can be observed that the results of PB3 method are very close to the exact solution with a small deviation. The average deviations from the optimality for these algorithms are depicted in Figure 6, the solution obtained from PB3 algorithm for samples (1) to (4) are the same as the exact solution, as their deviations are zero. Furthermore, in comparison with the PB2 and PB3 methods, PB3 method resulted in smaller deviation values as the size of the problem increases. The optimality deviation in problem instance #14 is 760.9 (0.0024%) for PB3 algorithm while it is 1856631.6 (5.773%) and 1354665.6 (4.213%) units for PB1 and PB2 methods, respectively.Figure 6
The average optimality deviation of the solutions obtained by each solution method within the small and medium size problem instances.To demonstrate the difference between the output of the mentioned algorithms and the exact solution obtained from CPLEX, the solution of the problem instance (1) is reported in Table7. According to this table, it is observed that all the three methods present the same number and locations for facilities. However, the difference is related to the quantity flow between facilities. The obtained quantity flow among facilities in PB3 methods and the exact solution is the same.Table 7
The location and quantity flow between facilities in problem instance #1.
Variables PB1 PB2 PB3 CPLEX X1,1 133.96 133.96 133.96 133.96 X3,2 126.2 126.2 126.2 126.2 X3,3 104.18 104.18 104.18 104.18 X3,4 128.34 128.34 128.34 128.34 X1,5 135.56 135.56 135.56 135.56 X3,6 127.21 127.21 127.21 127.21 X3,7 118.03 118.03 118.03 118.03 X1,8 86.95 126.21 126.21 126.21 X1,9 0 96.95 0 0 X3,8 39.26 0 0 0 X3,9 120.49 33.55 120.49 120.49 X3,10 116.73 116.73 116.73 116.73 T 1,1 236.4 236.4 236.4 236.4 Y 1,1 15.32 15.32 15.32 15.32 Y 2,1 11.22 11.22 11.22 11.22 Y 3,1 37.27 37.27 37.27 37.27 Y 4,1 42.78 42.78 42.78 42.78 Y 5,1 31.34 31.34 31.34 31.34 Y 6,1 25.15 25.15 25.15 25.15 Y 7,1 32.58 32.58 32.58 32.58 Y 8,1 48.33 48.33 48.33 48.33 Y9,1 39.11 39.11 39.11 39.11 Y 10,1 12.82 12.82 12.82 12.82 S 1,1 59.54 59.54 59.54 59.54


### 6.4. Convergence and Speed of Algorithms

To show the convergence rate of the algorithms and their quality, one of the problem instances, namely, instance #12, is used to demonstrate the performance of the proposed algorithms in Figure7. The optimality gap is shown for each of the iterations and it is clearly observable that PB3 method is able to further search the solution space, which raise the chance of obtaining the optimal solution.Figure 7
The gap between the optimal objective value and the output of the three algorithms over the iterations for a problem instance.The computational time of the algorithms is another interesting factor to point out specially in NP-hard problems. The average computational time of our algorithms for each problem instance is given in Figure8. According to this figure it is observed that the computational time of the exact solution increases exponentially as the size of the problem increases. The computational time of the exact solution has become more than that of PB1, PB2, and PB3 methods from the problem instance #3 onwards. This time increase shows that the exact solution is inefficient in solving large-scale problems and even it may not achieve the optimal solution. However, the computational times of PB1, PB2, and PB3 algorithms are close to each other and they are less than the solver’s except for some very small size problem instances. Generally, the study revealed that the proposed solution method has obtained results very close to those of the exact solution than the other methods. However, its computational time at large sizes is less than that of the exact solution.Figure 8
The average computational time of methods for each sample problem.


### 6.5. Algorithms Performance for Larger Scale Problems

To evaluate the performance of the proposed solution method for large-scale instances, 7 sample problems were solved using PB1, PB2, and PB3 algorithms each of them five times. As the solver could not find a feasible solution, the performance comparison on the solution quality is made using RPD measure which is defined below:(20)%RPD=Solmethod-BestSolBestSol×100, where Solmethod is the solution obtained by each algorithm and BestSol is the best solution obtained among all three solution methods. Table 8 presents the computational results in large-scale problems.Table 8
The average of the objective values for each problem instance (in large scale problems).
Problem instance # PB1 PB2 PB3 15 36913935 35518893 33938760 16 40865814 39330392 36800070 17 45761578 44397034 42238210 18 49926096 48933683 46045200 19 58336893 56463500 53420906 20 64568402 62410560 58308146 21 69765594 67512119 64137591Figures9-10 present a schematic comparison of the average CPU time and the average RPD value of three solution methods, respectively.Figure 9
The average CPU-Time of each solution method in each large-scale problems.Figure 10
The average RPD of each solution method in each large-scale problems.


### 6.6. Statistical Analysis

We have used Duncan’s multiple range test as a statistical technique to compare the average results of each problem instance to conclude about the performance of the solution methods. The steps of Duncan’s multiple range test are as follows:(1)
Create a tree diagram, whose first node includes the average ofK treatments. Then, sort them in an increasing form and call it p level.(2)
Create two branches from the first node. The first node includes the average of treatment 1 toK-1, and the second node includes the average of treatment 2 to K. Then, call it p-1 level.(3)
Proceed with step (2) pattern to the point that there are only two pairs of each node.(4)
Calculate the range ofR for each branch of p level in the tree chart using the following equation:(21)Rp=Serα,p,dfen,where Se is the standard error in the ANOVA, n is the number of observations in each treatment, and the values of rα,p,dfe are obtained from Duncan’s multiple range test table for the corresponding α value. rα,p,dfe depends on the significance level α, the number of treatments p at the node, and the degrees of freedom of the error, dfe, for the ANOVA. (1)
Start from the first node and compare theR value with the least significant level of Rp Duncan.
(a)
IfR>Rp, it can be concluded that the two averages are significantly different at the beginning and end of p level node. Thus, apply this test to the next nodes, which is shown by Yes.(b)
IfR<Rp, it can be concluded that there is no significant difference between the averages of p level node. Thus, do not apply this test to the next nodes, which is shown by No.In our test, the aforementioned solution methods were compared via 21 sample problems and the results are reported in Table9. For example, in the problem instance #12, it can be seen that there is a significant difference between the averages obtained from application of PB3 and PB2 methods while there is no significant difference between the averages obtained from application of PB2 and PB1 methods.Table 9
The results of the Duncan test among each solution method.
Problem instance # PB1 v.s. PB2 PB1 v.s. PB3 PB2 v.s. PB3 1 NO NO NO 2 NO NO NO 3 NO NO NO 4 NO NO NO 5 NO YES YES 6 NO YES YES 7 NO YES NO 8 NO NO YES 9 NO NO NO 10 NO YES YES 11 NO YES YES 12 NO YES YES 13 NO YES YES 14 NO YES YES 15 NO YES YES 16 NO YES YES 17 YES YES YES 18 YES YES YES 19 YES YES YES 20 YES YES YES 21 YES YES YES


## 6.1. Problem Instances

In order to assess the performance of the proposed discrete LCA in terms of the objective-function value and CPU time, several numerical experiments with different problem sizes were implemented, and the obtained results are reported in this section. To this end, 21 sample problems with 3 different levels, which are small, medium, and large-scale problems, were generated by different combinations of the parameter values. The levels of the generated sample problems are shown in Table3, and the ranges of the parameters are presented in Table 4. Furthermore, all parameters of the sample problems were randomly generated based on uniform distributions in prespecified intervals.Table 3
Size and scale of sample problems.
Scale of the problem Problem no. # the potential locations for plants # the potential locations for collection centers # the potential locations for disposal centers # the potential locations for customer zones Small 1 3 2 2 10 2 4 3 3 15 3 5 4 3 20 4 7 5 4 25 5 8 6 5 30 6 10 8 7 35 7 12 8 7 40  Medium 8 15 9 8 45 9 18 10 9 50 10 21 12 10 55 11 25 15 12 60 12 30 20 15 70 13 35 25 18 80 14 40 32 20 90  Large 15 45 37 24 100 16 50 42 29 120 17 60 48 33 140 18 70 53 38 160 19 80 57 42 180 20 90 60 45 200 21 100 65 50 220Table 4
Pre-specified intervals to generate parameters intervals based on uniform distributions.
Parameter Range Parameter Range Ei U ( 1000000,1200000 ) TSCjk U(20,30) Fk U ( 1000000,1200000 ) TCPki U(20,30) Gl U ( 1000000,1200000 ) TCDkl U(20,30) D e m j U(100,150) C a p P i U(800,1200) rj U ( 10,50 ) CapCk U(200,400) β U(0.6,0.8) CapDl U(100,300) T P S i j U(20,30)As the acquired results obtained from the proposed algorithm are sensitive to their initial parameters, any small changes can affect the accuracy of the best obtained solution. Therefore, the Taguchi tuning method is used for the parameters to find the best solutions. In this method, first, the appropriate factors (initial parameters) are determined, the level of each factor is selected, and then the design of experiments for this control factor is specified. After specifying the experimental design, the proposed algorithm is used to find the best combination of factors. In our study, 4 levels are considered for each factor. In Table5, the number of teams (l), number of seasons (S), and scaling constants of strength and weakness components (ψ1 and ψ2) are given.Table 5
Proposed parameter levels.
Parameter Level 1 Level 2 Level 3 Level 4 l 100 150 200 25 S 20 26 32 40 ψ1 0 2 4 6 ψ2 0 2 4 6


## 6.2. Objective Function Average

With respect to the proposed algorithm, the experimental design is employed according to the number of factors and number of levels. To this aim, each problem instance was solved five times to form the replications. Average results are reported as the final value. The best values of the proposed parameters for discrete LCA according to the mean normalized objective values are 150, 20, 6, and 2 for the number of teams, number of seasons, and constant coefficients to scale the contribution of the strength and weakness components, respectively. Moreover, Figure5 presents the normalized average of means and the S/N (Signal-to-noise) ratio graph for the experimental design of the algorithm. Based on the average of the means given in the graph, the algorithms result in a more efficient response if the parameter is placed in the lower level. Furthermore, in the S/N ratio graph, the algorithms result in a more efficient response if the parameter is placed in the higher level. According to Table 5 and Figure 5 the best performance of our algorithm is observed when the number of teams (l) and the weakness scaling coefficient (ψ2) are in level 2; the number of seasons (S) is placed at level 1, and the strength scaling coefficient (ψ4) is placed at level 4.Figure 5
Means of means and the S/N ratio plot for discrete LCA.


## 6.3. Solutions Structure and Quality

After tuning the parameters of the proposed algorithm, PB3 results are compared with those of the CPLEX outputs which were run on GAMS software (exact solution), PB2, and PB1 methods. For this purpose, five problem instances with random data given in Table4 for each problem size given in Table 3 were generated, and they were solved by each method. The average of five problems were selected as the benchmark for each sample problem. In Table 6, the average results of each sample problem instance is provided and compared with the solver’s (exact) solution.Table 6
The average of the objective values for each sample problem instance.
Scale of the problem Problem instance CPLEX PB1 PB2 PB3 Small 1 5216131.55 5216258.95 5216878.01 5216131.55 2 6537367.23 6537605.34 6538131.72 6537367.23 3 7709177.89 7709634.54 7710523.47 7709177.89 4 9378103.61 9378752.72 9379025.74 9378103.61 5 11611170.7 11612119.3 11612637.7 11611172.3 6 12217294.9 12219232.9 12218828.6 12217296.1 7 14546224.7 14547451.7 14547694.5 14546227.6  Medium 8 15642700.4 15645397.6 15644410.2 15642705.3 9 16394986.2 16411886.8 16408650.5 16395010.8 10 18192348.5 18209196.2 18199012.0 18192397.3 11 19613086.9 19800951.6 19628863.7 19613153.8 12 22984000.0 23693969.4 23485221.2 22984540.7 13 25164851.6 26212668.5 26139382.5 25165442.0 14 32156898.7 34013530.3 33511564.3 32157659.6According to the results shown in Table6, it can be observed that the results of PB3 method are very close to the exact solution with a small deviation. The average deviations from the optimality for these algorithms are depicted in Figure 6, the solution obtained from PB3 algorithm for samples (1) to (4) are the same as the exact solution, as their deviations are zero. Furthermore, in comparison with the PB2 and PB3 methods, PB3 method resulted in smaller deviation values as the size of the problem increases. The optimality deviation in problem instance #14 is 760.9 (0.0024%) for PB3 algorithm while it is 1856631.6 (5.773%) and 1354665.6 (4.213%) units for PB1 and PB2 methods, respectively.Figure 6
The average optimality deviation of the solutions obtained by each solution method within the small and medium size problem instances.To demonstrate the difference between the output of the mentioned algorithms and the exact solution obtained from CPLEX, the solution of the problem instance (1) is reported in Table7. According to this table, it is observed that all the three methods present the same number and locations for facilities. However, the difference is related to the quantity flow between facilities. The obtained quantity flow among facilities in PB3 methods and the exact solution is the same.Table 7
The location and quantity flow between facilities in problem instance #1.
Variables PB1 PB2 PB3 CPLEX X1,1 133.96 133.96 133.96 133.96 X3,2 126.2 126.2 126.2 126.2 X3,3 104.18 104.18 104.18 104.18 X3,4 128.34 128.34 128.34 128.34 X1,5 135.56 135.56 135.56 135.56 X3,6 127.21 127.21 127.21 127.21 X3,7 118.03 118.03 118.03 118.03 X1,8 86.95 126.21 126.21 126.21 X1,9 0 96.95 0 0 X3,8 39.26 0 0 0 X3,9 120.49 33.55 120.49 120.49 X3,10 116.73 116.73 116.73 116.73 T 1,1 236.4 236.4 236.4 236.4 Y 1,1 15.32 15.32 15.32 15.32 Y 2,1 11.22 11.22 11.22 11.22 Y 3,1 37.27 37.27 37.27 37.27 Y 4,1 42.78 42.78 42.78 42.78 Y 5,1 31.34 31.34 31.34 31.34 Y 6,1 25.15 25.15 25.15 25.15 Y 7,1 32.58 32.58 32.58 32.58 Y 8,1 48.33 48.33 48.33 48.33 Y9,1 39.11 39.11 39.11 39.11 Y 10,1 12.82 12.82 12.82 12.82 S 1,1 59.54 59.54 59.54 59.54


## 6.4. Convergence and Speed of Algorithms

To show the convergence rate of the algorithms and their quality, one of the problem instances, namely, instance #12, is used to demonstrate the performance of the proposed algorithms in Figure7. The optimality gap is shown for each of the iterations and it is clearly observable that PB3 method is able to further search the solution space, which raise the chance of obtaining the optimal solution.Figure 7
The gap between the optimal objective value and the output of the three algorithms over the iterations for a problem instance.The computational time of the algorithms is another interesting factor to point out specially in NP-hard problems. The average computational time of our algorithms for each problem instance is given in Figure8. According to this figure it is observed that the computational time of the exact solution increases exponentially as the size of the problem increases. The computational time of the exact solution has become more than that of PB1, PB2, and PB3 methods from the problem instance #3 onwards. This time increase shows that the exact solution is inefficient in solving large-scale problems and even it may not achieve the optimal solution. However, the computational times of PB1, PB2, and PB3 algorithms are close to each other and they are less than the solver’s except for some very small size problem instances. Generally, the study revealed that the proposed solution method has obtained results very close to those of the exact solution than the other methods. However, its computational time at large sizes is less than that of the exact solution.Figure 8
The average computational time of methods for each sample problem.


## 6.5. Algorithms Performance for Larger Scale Problems

To evaluate the performance of the proposed solution method for large-scale instances, 7 sample problems were solved using PB1, PB2, and PB3 algorithms each of them five times. As the solver could not find a feasible solution, the performance comparison on the solution quality is made using RPD measure which is defined below:(20)%RPD=Solmethod-BestSolBestSol×100, where Solmethod is the solution obtained by each algorithm and BestSol is the best solution obtained among all three solution methods. Table 8 presents the computational results in large-scale problems.Table 8
The average of the objective values for each problem instance (in large scale problems).
Problem instance # PB1 PB2 PB3 15 36913935 35518893 33938760 16 40865814 39330392 36800070 17 45761578 44397034 42238210 18 49926096 48933683 46045200 19 58336893 56463500 53420906 20 64568402 62410560 58308146 21 69765594 67512119 64137591Figures9-10 present a schematic comparison of the average CPU time and the average RPD value of three solution methods, respectively.Figure 9
The average CPU-Time of each solution method in each large-scale problems.Figure 10
The average RPD of each solution method in each large-scale problems.


## 6.6. Statistical Analysis

We have used Duncan’s multiple range test as a statistical technique to compare the average results of each problem instance to conclude about the performance of the solution methods. The steps of Duncan’s multiple range test are as follows:(1)
Create a tree diagram, whose first node includes the average ofK treatments. Then, sort them in an increasing form and call it p level.(2)
Create two branches from the first node. The first node includes the average of treatment 1 toK-1, and the second node includes the average of treatment 2 to K. Then, call it p-1 level.(3)
Proceed with step (2) pattern to the point that there are only two pairs of each node.(4)
Calculate the range ofR for each branch of p level in the tree chart using the following equation:(21)Rp=Serα,p,dfen,where Se is the standard error in the ANOVA, n is the number of observations in each treatment, and the values of rα,p,dfe are obtained from Duncan’s multiple range test table for the corresponding α value. rα,p,dfe depends on the significance level α, the number of treatments p at the node, and the degrees of freedom of the error, dfe, for the ANOVA. (1)
Start from the first node and compare theR value with the least significant level of Rp Duncan.
(a)
IfR>Rp, it can be concluded that the two averages are significantly different at the beginning and end of p level node. Thus, apply this test to the next nodes, which is shown by Yes.(b)
IfR<Rp, it can be concluded that there is no significant difference between the averages of p level node. Thus, do not apply this test to the next nodes, which is shown by No.In our test, the aforementioned solution methods were compared via 21 sample problems and the results are reported in Table9. For example, in the problem instance #12, it can be seen that there is a significant difference between the averages obtained from application of PB3 and PB2 methods while there is no significant difference between the averages obtained from application of PB2 and PB1 methods.Table 9
The results of the Duncan test among each solution method.
Problem instance # PB1 v.s. PB2 PB1 v.s. PB3 PB2 v.s. PB3 1 NO NO NO 2 NO NO NO 3 NO NO NO 4 NO NO NO 5 NO YES YES 6 NO YES YES 7 NO YES NO 8 NO NO YES 9 NO NO NO 10 NO YES YES 11 NO YES YES 12 NO YES YES 13 NO YES YES 14 NO YES YES 15 NO YES YES 16 NO YES YES 17 YES YES YES 18 YES YES YES 19 YES YES YES 20 YES YES YES 21 YES YES YES


## 7. Conclusion and Suggestions for Future Studies

In this paper, a CLSC network including levels of the plants to customer zones, customer zones to collection centers, collection centers to plants, and collection centers to disposal centers was considered. We proposed a mathematical programming model and then developed a new population-based algorithm called the LCA to solve the model.The proposed solution algorithm was designed by modifying an existing LCA in the literature which is for continuous space. We defined new operators and used them to convert the continuous space into the discrete one and hence we named our algorithm discrete LCA. We also proposed a new priority-based encoding for solving the problems. To demonstrate the effectiveness of the proposed method, an experimental study was conducted to study the small, medium, and large-scale problem instances. The numerical study showed that the outputs of our solution method is very close to the exact optimal solution which are obtained from the solver while it is significantly faster for large size problem instances. Furthermore, from the aspect of objective value, our proposed encoding method outperforms the two other existing methods in the literature which we referred to as PB1 and PB2 within the text whereas there are no significant differences between their computational time.This work can be extended either with augmenting the model assumption or from the solution approach. For instance, considering a multiobjective model which addresses the sustainability concerns in supply chain network design, either from the customer viewpoint or from the suppliers, would be an interesting problem to investigate in this framework. In addition, considering uncertainty in the demands and the returned-products can result in a more realistic and challenging problem to deal with.



---

*Source: 8163927-2018-05-17.xml*